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1.  INTRODUCTION 


1.1  Background 

In  the  search  for  UXO  and  for  discrimination  between  UXO  and  non-UXO  metallic 
fragments  (clutter)  it  is  necessary  to  accurately  determine  the  parameters  that  characterize  a 
metallic  object  in  the  ground.  A  search  system  is  needed  that  not  only  detects  the  object 
itself  but  can  also  determine  its  size,  shape,  orientation,  shell  thickness  and  metal  content 
(ferrous  or  non-ferrous,  mixed  metals).  These  properties  of  a  buried  metallic  object  are 
referred  to  here  as  the  object  parameters. 

The  search  for  UXO  is  a  two-step  process.  The  object  must  first  be  detected  and  its 
location  determined  then  the  parameters  of  the  object  must  be  defined.  The  first  step  is 
now  accomplished  with  a  variety  of  magnetometer  and  active  electromagnetic  (AEM) 
systems.  The  AEM  systems  operate  in  the  transient  or  frequency  domain  mode  and  at 
present  use  a  single  transmitter  and  up  to  three  receivers.  The  characterization  phase  can 
also  be  described  as  a  two-step  process.  A  variety  of  incident  fields  are  used  to  induce 
magnetization  and  current  flow  in  different  directions  within  the  object.  The  magnetic 
dipole  moments  induced  in  the  body,  normalized  by  the  inducing  field  are  known  as  the 
polarizabilities  of  the  object.  The  secondary  fields  related  to  these  induced  polarizabilities, 
as  a  function  of  frequency  or  time,  are  the  measured  quantities.  Interpretation  could  stop 
with  the  determination  of  the  object  location  and  the  three  orthogonal  principal 
polarizabilities.  These  polarizabilities  and  their  variation  with  either  time  or  frequency  are 
the  only  fundamental  object  parameters  that  can  be  recovered  from  the  inductive  excitation 
of  a  finite  body  in  the  ground  if  a  dipolar  representation  is  assumed.  Presumably  a  catalog 
of  polarizabilities  could  be  constructed  from  the  forward  model  polarizabilities  for  a  wide 
range  of  potential  targets  and  then  a  look-up  table  inversion  scheme  could  be  used  to 
identify  the  actual  dimensions  and  physical  properties  of  a  target. 

We  have  found  that  a  powerful  second  step  in  the  characterization  process  is  to 
directly  interpret  the  polarizabilities  in  terms  of  the  principal  axes  of  an  equivalent 
spheroid.  This  step  yields  the  size  and  aspect  ratio  of  the  spheroid  that  is  physically 
equivalent  to  the  target.  Finally,  the  wideband  target  response  permits  the  estimation  of 
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both  the  permeability-conductivity  ratio  and  the  conductivity-permeability  product  so  that 
each  can  be  estimated  independently.  It  is  unlikely  that  any  better  resolution  of  the  shape 
of  an  intact  UXO  than  its  equivalent  spheroid  can  be  recovered  using  practical  system  data. 
A  catalog  of  equivalent  spheroidal  bodies  representing  most  UXOs  would  thus  be  relatively 
easy  to  construct  for  final  target  identification.  This  more  detailed  second  step  requires  the 
broadband  spectral  or  transient  response  using  frequencies  low  enough  to  identify  the  quasi 
dc  magnetization  response  and  high  enough  to  identify  the  purely  electromagnetic  (EM), 
inductive  limit,  response  which  depends  only  on  the  size  of  the  object. 

The  concept  of  characterizing  UXO  by  equivalent  spheroid  parameters  is  key  to 
distinguishing  intact  UXO  from  non-UXO  metal  scrap.  Any  UXO  is  expected  to  retain  its 
fundamental  shape  (size,  aspect  ratio  and  symmetry  about  its  long  axis)  with  perhaps  minor 
distortion  caused  by  impact.  Metal  scrap  will  have  distinct  polarizability  signatures  that 
cannot  mimic  these  of  elongated  symmetric  bodies.  Roughly  flat  sheets  will  have  dipolar 
responses  approaching  these  of  a  highly  flattened  oblate  spheroid  (close  to  a  loop 
response),  twisted  sheets  a  principal  moment  orthogonal  to  some  equivalent  plane  through 
the  sheet  with  small  and  highly  irregular  minor  axis  polarizabilities  etc.  These 
distinguishing  polarizabilities,  coupled  with  the  size  estimates  and  spatial  sampling  of  the 
multiple  receiver  array  are  more  than  enough  to  separate  small  scrap  from  deeper  targets  of 
concern. 

An  optimal  EM  system  can  extract  from  the  measurements  the  best  possible 
estimates  of  the  location,  orientation,  size,  shape  and  metal  content  of  a  buried  metallic 
object  in  the  presence  of  the  interfering  response  of  the  ground  and  non-UXO  metallic 
objects.  Discrimination  can  be  achieved  partly  by  selective  filtering  of  the  responses 
inherent  in  the  system  configuration  and  design  and  partly  through  post  acquisition 
processing.  The  target  parameters  are  obtained  by  a  statistical  inversion  of  the 
measurements  to  establish  the  principal  electromagnetic  moments  of  the  detected  object 
and  their  variation  with  time  or  frequency. 
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1.2  Objectives 


In  summary  the  technical  objectives  of  this  project  were: 

1)  To  develop  and  demonstrate  a  methodology  for  the  quantitative  evaluation  of  existing 
AEM  systems  and  for  the  design  of  new  systems. 

2)  To  implement  a  new  methodology  for  optimizing  an  AEM  system  for  detecting  and 
classifying  UXO  of  a  given  class  in  a  specified  geologic  setting  and  in  a  given  noise 
environment. 

3)  To  design  and  build  a  prototype  of  an  active  EM  system  for  detecting  and  characterizing 
a  metallic  object  in  the  ground. 

1.3  Technical  Approach 

The  design  methodology  is  based  on  the  use  of  simulators,  numerical  models  of  the 
electromagnetic  response  of  an  arbitrary  target  in  the  ground  to  an  arbitrary  configuration 
of  transmitters  and  receivers.  The  variables  in  the  simulator  are  the  parameters  of  the 
targets,  the  parameters  of  the  AEM  system  and  the  noise  (ambient,  motion,  instrumental, 
and  geologic).  The  target  parameters  are  the  location,  orientation,  depth,  size,  shape,  metal 
content  and  type.  A  schematic  of  the  simulator  is  shown  in  Figure  1.3.1. 


EVALUATION 


Figure  1.3.1:  Schematic  diagram  of  the  simulator. 
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OPTIMIZATION 


l 


OPTIMUM  SYSTEM 
PARAMETERS 


Figure  1.3.2:  Schematic  diagram  of  the  optimization  process. 


VERIFICATION 


Figure  1.3.3:  Schematic  diagram  of  the  verification  process. 


The  simulator  can  be  used  to  evaluate  designs  by  simulating  the  response  of  a  given 
system,  with  a  given  noise  level,  to  a  particular  target.  Various  system  configurations  can 
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be  quantitatively  compared  through  their  respective  signal  to  noise  ratios  over  the  same 
target. 

This  forward  modeling  approach  is  useful  for  evaluating  the  relative  response  of 
different  targets  for  a  given  system,  for  example  in  analyzing  the  role  of  loop  size  in 
discriminating  between  small,  shallow,  targets  and  deeper,  larger,  ones.  It  is  also  useful  for 
evaluating  the  ground  response  and  modeling  the  spectral  or  transient  response  of  various 
targets. 

The  system  parameters  that  are  variables  in  the  design  are:  a)  the  geometric 
configuration  of  the  transmitter(s)  and  receivers(s)  and  the  way  in  which  they  are  mounted 
for  given  target  objectives  (the  platform),  b)  the  spatial  positioning  of  the  system  (profile, 
grid,  single  site  stand-off,  etc.),  c)  the  transmitter  power  and  waveform,  d)  the  system  and 
ambient  noise,  e)  the  receiver  bandwidth  and  dynamic  range,  f)  the  signal  averaging  time  (a 
function  of  survey  speed).  The  geologic  variables  (geologic  noise)  are  the  values  and 
variability  of  the  ground  conductivity  and  permeability.  A  more  detailed  schematic  of  all 
the  variables  that  contribute  to  the  data  acquired  with  a  general  AEM  system  is  shown  in 
Figure  1.3.4. 


Figure  1.3.4:  Generic  AEM  System. 
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In  inverse  mode  the  simulator  is  used  to  determine  target  properties  from 
measurements  with  a  given  transmitter-receiver  configuration.  At  any  given  time  in  the 
transient  response  (or  at  a  given  frequency  in  the  frequency  domain)  the  response  data  are 
inverted  to  yield  the  location  (x,  y,  z)  of  the  target,  its  attitude  and  its  principal 
polarizabilities  (yielding  an  apparent  aspect  ratio).  In  this  process  noise  estimates  (or 
measurements)  are  interpreted  to  yield  error  estimates  of  the  object  location,  attitude  and 
polarizabilities.  This  inversion,  at  a  succession  of  times  in  the  transient,  or  at  many 
frequencies,  yields  the  object  polarizabilities  as  a  function  of  time  or  frequency,  which  can 
in  turn  yield  the  size,  true  aspect  ratio  and  estimates  of  the  conductivity  and  permeability  of 
the  target.  The  accuracy  of  these  property  estimates  depends  on  the  time  or  frequency 
window  over  which  the  polarizability  measurements,  and  their  accuracies,  are  known.  This 
particular  inversion  process  can  be  used  on  artificial  test  data  to  determine  the  bandwidth  of 
the  system  that  is  needed  to  provide  the  desired  accuracy  in  the  physical  property 
measurement. 

Finally,  the  simulator  was  used  to  determine  the  optimum  configuration  for 
detecting  and  characterizing  the  target.  In  this  step  certain  transmitter  configurations  are 
selected  and  the  position  and  orientation  of  n  receivers  (with  given  signal-to-noise  ratios) 
are  varied  until  targets  of  various  shape  and  depth  are  characterized.  The  array  is  judged  as 
optimum  when  the  errors  in  target  parameter  estimation  are  minimized.  An  initial  analysis 
using  this  approach  for  targets  of  specific  polarizability  led  to  the  conclusion  that  a  system 
employing  three  orthogonal  transmitters  and  a  minimum  of  five  vertical  receivers  deployed 
within  the  footprint  of  the  horizontal  loop  transmitter  would  be  optimal.  In  this  analysis, 
the  criterion  for  optimization  was  that  the  sum  of  the  squared  error  in  the  individual 
principal  polarizability  be  a  minimum.  This  left  open  the  possibility  that  certain 
configurations  might  be  used  to  estimate  one  polarizability  very  well  while  yielding  poor 
estimates  of  the  other  two. 

Since  the  accurate  determination  of  the  polarizabilities  is  needed  to  determine  the 
true  aspect  ratio  and  physical  properties  of  the  target,  the  inversion  algorithm  was  extended 
to  determine  the  number  of  receivers  and  their  orientation  needed  to  resolve  polarizabilities 
of  any  arbitrary  target  located  beneath  the  transmitter-receiver  array.  This  analysis  showed 
that  more  than  five  receivers  would  be  advantageous  but  that  using  multiple  receiver 
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orientations  had  little  practical  advantage  over  slightly  larger  arrays  of  only  vertical 
receivers. 

The  results  of  all  these  quantitative  analyses  led  to  the  preliminary  design  of  a 
configuration  using  three  orthogonal  loop  transmitters  and  eight  vertical  receivers  deployed 
in  an  asymmetric  pattern  in  the  plane  of  the  horizontal  loop.  This  is  shown  schematically 
in  Figure  1.3.5  and  a  photo  of  the  bench  prototype,  showing  one  of  the  loop  transmitters,  is 
shown  in  Figure  1.3.6. 


target 


Figure  1.3.5:  Schematics  of  the  bench  test  acquisition  setup. 


Figure  1.3.6:  Photo  of  the  prototype  system. 
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2.  INVERSION  FOR  TARGET  PARAMETERS  AND  THE  ROLE  OF  SYSTEM 

CONFIGURATION 


2.1  Estimating  Equivalent  Dipole  Polarizabilities  for  the  Inductive  Response  of 
Isolated  Conductive  Bodies 

2.1.1  Introduction 

Any  set  of  currents  can  be  characterized  in  terms  of  a  set  of  their  multipole 
moments.  The  associated  magnetic  fields  can  be  represented  as  a  sum  of  corresponding 
multipole  terms  away  from  the  currents  (e.g.,  Jackson,  1975,  p.746).  For  a  magnetic 
multipole  term  of  order  n,  magnetic  field  strengths  fall  off  as  1/r  n+2,  in  resistive  media. 
Dipole  terms  are  the  lowest  order  magnetic  multipole  terms.  At  distances  much  greater 
than  the  scale  of  an  object,  dipole  terms  become  a  very  good  approximation  to  the 
magnetic  fields  arising  from  currents  induced  in  the  object. 

In  the  vicinity  of  a  conductive  body,  the  primary  magnetic  field  imposed  by  an 
external  source  current  may  be  approximated  by  the  primary  magnetic  field  at  the  object 
center  ro,  B(p)(ro,t).  Assuming  a  common  time  variation  g(t)  for  all  primary  magnetic 
field  components  at  the  object  center,  we  define  a  primary  field  magnitude  vector  as 
B(0)  =  B<p)(r0,t)/g(t).  We  choose  the  normalization  of  g(t)  so  that  g(to)  =  1  at  some  chosen 
time  to,  for  example,  for  a  step  function  turn-off  primary  field  we  choose  the  scale  of  g(t) 
so  that  g(t)  =  1  for  t<0.  In  practice,  it  is  common  to  assume  that  the  medium  surrounding 
the  object  is  sufficiently  resistive  that  the  magnetic  fields  due  to  currents  induced  in  the 
surrounding  medium  are  negligible  at  the  body,  so  that  g(t)  is  simply  the  transmitter 
current  waveform.  Neglecting  primary  field  gradients,  the  secondary  magnetic  fields 

B(s)(r,t)  =  B(r,t)  -  B(p)(r,t)  (2.1.1) 

due  to  currents  induced  in  a  conductive  body  can  be  written  as  linear  combinations  of  the 
magnetic  fields  that  would  be  induced  by  primary  fields  of  strength  g(t)  in  the  x,  y,  or 

£  direction  at  the  objects  center,  Bx(s\r,t),  By<s)(r,t),  or  Bzls,(r,t)  respectively; 

B(s)(r,t)  =  Bx(0)  Bx(s)(r,t)  +  By(0)By(s)(r,t)  +  Bz(0)  Bz(s)(r,t)  ,  (2. 1 .2) 

where  Bx(0),  By(0),  Bz(0)  are  the  x,  y,  and  z  components  of  Bl0).  Being  secondary  magnetic 
fields  induced  by  a  source  magnetic  field  of  strength  g(t)  (which  is  dimensionless) 
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Bx(s)(r,t),  By<s,(r,t),  and  Bz(s)(r,t)  are  similarly  dimensionless.  Forming  a  matrix  with  the 
three  vectors  Bx(s|(r,t),  By(s)(r,t),  and  Bz<s)( r,t)  as  its  columns,  the  right  hand  side  of 
Equation  2. 1 .2  is  a  matrix  vector  product; 

B(s)(r,t)  =  [Bx(s)(r,t),  By(s)(r,t),  Bz(s)(r,t)]  •  B(0)  ,  (2. 1 .3) 

where  B<s,(r,t),  Bx(s)(r,t),  By(s)(t\t),  Bz(s)(r,t)  and  B<0)  considered  as  column  vectors,  the  dot 
effect  matrix  multiplication. 

At  distances  where  non-dipole  secondary  magnetic  fields  are  small,  the  secondary 
magnetic  fields  induced  by  the  primary  magnetic  field  in  the  x  direction  can  be  broken 
into  contributions  by  dipole  components  in  the  x,  y,  and  z  directions; 

Bx(s)(r,t)  =  mxx(t)  Bx<d)(r)  +  myx(t)By(d)(r)  +  mzx(t)  Bz(d)(r),  (2.1 ,4a) 

where  Bx(d)(r),  By(d)(r),  and  B  z(d)(r)  are  the  magnetic  fields  of  a  unit  magnetic  dipole  in 
the  x,  y,  and  £  directions  respectively,  placed  at  the  body  center,  and  have  units  of 
Tesla/ Amp-m  .  Quantities  mxx(t),  myx(t),  and  mzx(t)  are  the  current  dipole  moments  in 
these  directions,  for  a  unit  primary  (inducing)  magnetic  field  in  the  x  direction  at  the 
object  center,  and  have  units  of  Amp-m2/Tesla.  Similarly, 

By(s)(r,t)  =  mxy(t)  Bx(d)(r)  +  myy(t)By(d)(r)  +  mzy(t)  Bz(d)(r),  (2. 1 ,4b) 

Bz(s)(r,t)  =  mxz(t)  Bx<d)(r)  +  myz(t)By(d)(r)  +  mzz(t)  Bz(d)(r),  (2.1 .4c) 

with  mxy(t),  myy(t),  mzy(t)  and  mxz(t),  myz(t),  mzz(t)  the  corresponding  moments  for 
primary  magnetic  fields  in  the  y  and  £  directions.  Assuming  that  the  surrounding 
medium  is  sufficiently  resistive  that  tertiary  currents  induced  in  the  surrounding  medium 
by  the  magnetic  fields  due  to  currents  in  the  body  can  be  neglected,  the  effective 
magnetic  dipole  moments  correspond  to  the  actual  moments  of  the  currents  circulating  in 
the  body.  We  make  this  assumption,  and  henceforth  refer  to  them  simply  as  the  dipole 
moments.  Equations  (2.1.4)  can  be  written  in  matrix  form  as 


Bx(s)(r),  B  (s)(r),  Bz(s)(r)l  =  [Bx(d)(r),  By(d)(r),  Bz(d)(r)' 


mxxmxymxz 


myxmyymyz 


mzxmzym2Z 


(2.1.5) 


where  the  explicit  time  dependence  of  the  matrix  of  dipole  moments  has  been  omitted. 
The  matrix  of  dipole  moments  M  is  symmetric  (Landau  and  Lifshitz,  1960,  pi  92). 
Substituting  equation  (2.1.5)  into  equation  (2.1.3)  gives 
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B(s)(r,t)  = 


Bx(d)(r),  B  (d)(r),  Bz(d)(r) 


M(t)  B1 


(0) 


(2.1.6) 


In  time  domain  applications,  M  is  real,  in  addition  to  being  symmetric,  so  can  be 
diagonalized  by  an  orthogonal  matrix  U(t); 

L(t)  =  UT(t)  M(t)  U(t)  (2.1 ,7a) 

where  T  denotes  transpose,  L(t)  is  diagonal,  with  elements  Ln(t),  L22(t),  L33  (t)  known  as 
the  principal  moments  of  M(t),  and  have  the  same  units  as  M  (Amp-m2/Tesla).  Equation 
(2.1.7a)  expresses  M(t)  in  coordinates  given  by  the  columns  of  U(t),  (u,),  known  as  the 
principal  directions  of  M(t).  For  bodies  with  an  axis  of  symmetry  w ,  w  is  one  principal 
direction  (e.g.,  Ui),  with  corresponding  principal  component  (e.g.,  Ln)  giving  the 
equivalent  dipole  moment  induced  in  the  w  direction  for  a  unit  primary  field  in  the  w 
direction  at  the  object  center.  The  other  two  principal  moments  correspond  to  equivalent 
dipole  moments  induced  in  directions  normal  to  w  for  unit  primary  fields  in  those 
directions.  Symmetry  of  the  object  implies  that  the  latter  two  moments  are  equal.  For  a 
symmetric  object,  rotating  into  coordinates  aligned  with  the  object's  symmetry  axis 
diagonalizes  M,  so  may  be  accomplished  by  a  rotation  matrix  U,  which  is  independent  of 
time. 


This  definition  of  the  equivalent  dipole  polarizability  matrix  M  is  consistent  with 
that  used  by  Pasion  and  Oldenburg  (2001),  and  differs  by  a  factor  of  po  from  that  used  by 
Baum  (1999). 

As  written,  equation  (2.1.6)  represents  the  magnetic  field  B(0)(r,t)  as  a  linear 
combination  of  dipole  fields.  Differentiating  it,  one  can  apply  the  same  methods  to 
modeling  measurements  of  d  B(s)  (r,  t)/dt,  with  dM(t)/dt  replacing  M(t).  In  this  case, 
dM(t)/dt  may  be  diagonalized  analogously  to  Equation  (2.1.7a); 

L’(t)  =  U’T(t)  M(t)  U’(t)  (2. 1 ,7b) 

Where  for  a  symmetric  object,  diagonalization  may  be  accomplished  with  the  same  time 
dependent  matrix  as  before;  U’(t)  =  U.  Strictly  speaking  dM/dt  and  L’(t)  represent 
equivalent  dipole  polarizability  (decay)  rates  rather  than  equivalent  dipole 
polarizabilities.  In  common  usage,  they  are  referred  to  simply  as  polarizabilities,  but  can 
be  distinguished  by  having  units  of  Amp-m  /s-Tesla  rather  than  Amp-m  /Tesla. 

For  a  given  time  dependence  of  source,  g(t),  equation  (2.1.6)  relates  secondary 
fields  at  any  time  to  an  equivalent  dipole  polarizability  M(t)  for  that  time,  so  M(t)  may  be 
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estimated  separately  for  each  time.  Consequently,  we  drop  the  explicit  time  dependence, 
and  assume  that  all  measurements  are  at  a  single  time  relative  to  the  starting  time  for  the 
primary  field  pulse. 


2,1.2  Estimating  Dipole  Polarizabilities  When  Object  Center  is  Known 

When  the  object  center  location  is  known,  B<0)  can  be  calculated  for  each  of  a  set 
of  sources  with  the  same  time  dependence  g(t).  For  the  i'th  measurement  of  a  set  of  n 
measurements,  letting  Bi(0)  be  the  primary  field  polarization  vector  at  the  object  center,  r, 
be  the  location  of  a  magnetic  field  measurement,  and  v(.  be  the  orientation  of  the  magnetic 


field  receiver  (e.g.,  coil),  equation  (2.1.6)  written  for  the  v:  component  at  r,  is 


v1TB(s)(r1)  =  v1T-  Bx d  (r;),  B®^),  Bz(d)(r,) 


M  B 


(0) 


(2.1.8) 


one  (scalar)  equation  constraining  the  six  unknown  dipole  polarizabilities  mxx,  niyy,  mzz, 
mxy  =  niyx,  myz  =  mzy,  and  mxz  =  mzx,  for  each  receiver  source  combination.  This  can  be 
rewritten  as 


di  f  i  xx  Hlxx  T  f  i  yy  ttlyy  f  i  zz  IWz  f  i  xy  Hfiy  +  f  i  yz  TXlyz  +  f  i  xz  ttlxz  5  (2. 1 .9) 

where  d,  =  viTB(s)(ri) ,  and  the  coefficients  f ,  xx ,  f ,  yy ,  •••  ,  can  be  found  by  multiplying 
out  the  vector  and  matrix  products  on  the  right  side  of  equation  (2.1.8),  substituting  mxy, 
myz,  and  mxz  for  niyx,  m7y,  and  mzx,  and  identifying  the  coefficients  of  mxx,  myy,  mzz,  mxy, 
myz,  and  mxz.  Equations  (2.1.9)  can  be  written  in  matrix  form  as 

d  =  F  m  ,  (2.1.10) 

where  m  =  (mxx  myy  ,  mzz ,  mxy  ,  myz  ,  mxz  )T ,  which  has  least  squares  solution 

m  =  (Ft  F)  1  Ft  d  .  (2.1.11) 

The  dipole  polarizability  moment  matrix  M  can  be  assembled  from  the  elements  of  m, 
using  the  symmetry  of  M. 


2.1.3  Estimating  Dipole  Polarizabilities  and  Object  Center  Location 

When  the  object  center  position  ro  is  unknown,  one  may  form  equation  (2.1.8) 
using  dipole  fields  Bx(d),  By(d),  BZW)  calculated  for  dipoles  centered  at  some  candidate 
object  center  position  r0,  and  primary  field  polarization  vectors  Bi<0)  at  the  candidate 
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object  center  position,  form  equation  (2.1.10),  and  calculate  the  least  squares  dipole 
polarizabilities  m  for  that  candidate  center  location,  m  (r0).  Its  squared  misfit  is 

X2  =|d-d(r0)|2/a2  ,  (2.1.12) 

where  a2  is  the  squared  measurement  uncertainty,  and 

d(r0)  =  F(FTF)_1FTd  (2.1.13) 

is  the  best  fitting  data  predicted  for  this  choice  of  ro.  Matrix  F  depends  on  ro  through  Bx 
(d),  By  (d),  Bz  (d),  and  B,  (0).  We  find  the  position  ro  giving  a  minimum  of  squared  misfit 
(2.1.12),  using  the  downhill  simplex  algorithm  (Press,  et  al.,  1986,  p.  289),  started  from 
four  candidate  center  locations,  r0®,  j=l,...,  4  at  the  comers  of  a  tetrahedron  with  edges 
one  quarter  of  the  length  of  the  maximum  separation  of  receiver  locations,  centered  one 
half  the  maximum  receiver  separation  below  the  receivers.  The  downhill  simplex 
algorithm  moves  the  corners  of  the  tetrahedron  systematically  expanding  or  contracting 
as  necessary  to  arrive  at  a  minimum  of  the  minimized  function,  (the  squared  misfit),  and 
ends  when  the  comers  have  converged  within  a  small  tolerance  of  each  other,  or  the 
function  values  at  the  four  comers  are  within  a  small  tolerance  of  each  other. 


2.1.4  Estimating  Dipole  Polarizability  Uncertainties 

For  data  with  small  measurement  errors,  the  uncertainty  in  the  resultant  dipole 
polarizabilities  and  equivalent  dipole  position  (object  center)  may  be  obtained  from 
analysis  of  a  linearized  inversion  for  dipole  polarizabilities  and  position.  Denoting  dipole 
components  at  the  i'th  receiver  in  the  receiver  direction  v,  by 


-  VB/Vi),  B'2(d)  =  VB/Vi),  B'i3  d  =  VB/Vi), 


(2.1.14) 


for  x,  y,  and  z  dipoles  respectively,  denoting  x,  y,  and  z  components  of  the  center  primary 
field  Bi<0)  by  B|,<0),  B2j<0),  and  B3/0),  and  numbering  the  elements  of  M  as  mkj,  for  k  =  1,  3, 
j=l,  3,  then  equation  (2.1.8)  can  be  written  explicitly  as 


di=Z2X 


(2.1.15) 


When  equivalent  dipole  position  r0  is  not  known  a  priori,  one  can  expand  equation  (2.15) 
in  a  Taylor  series  about  an  initial  value  r0(q),  such  as  the  result  of  the  downhill  simplex 
method  search  of  the  previous  section.  Letting  M(q)  be  the  corresponding  dipole 
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polarizability  matrix  fit  for  candidate  dipole  position  ro(q),  a  first  order  Taylor  expansion 
about  r0(q)  yields 

d,  =  +  mkj(<l)  [r0<<1+1)  -r0(<»]T  •  (2.1.16) 

k=l  j=l 

Collecting  coefficients  of  the  new  polarizability  estimates  mn(q+1),  m22(q+l) , into  a  row 
vector  a; <q  ),  and  coefficients  of  the  components  of  change  vector  Ar0  =  r0(q+l)  -  r0(C|) 
into  a  row  vector  gj(q)  equation  (2.1.16)  becomes 

di  =(ai(q),gi(q))  (m11,m22,m33,m12,m23,m13,Ax0,Ay0,Az0)T  ,  (2.1.17) 

where  superscript  (q+1)  has  been  omitted  from  the  various  nij/q 1  and  the  symmetry  of 
M  has  been  used  to  eliminate  ni2i(q+l),  m32(q+1),  and  m3i(q+1).  This  can  be  written  in  matrix 
form  as 

d  =  F  m  ,  (2.1.18) 

with  the  rows  of  F  and  vector  in  are  the  vectors  on  the  right  side  of  equation  (2.17),  and 
solved  for  min  the  same  manner  as  equations  (2.1.10)  and  (2.1.11).  A  new  estimated 
object  center  position  is  given  by 

i'„<,*')=r,<,>  +  Ar0  .  (2.1.19) 

Taylor  expanding  about  the  new  estimate  ro(q+1)  (equation  2.1.16,  with  q  incremented), 
the  process  is  repeated  until  the  change  magnitude  |Aro|  is  less  than  a  small  tolerance. 
The  variances  of  the  resultant  dipole  polarizabilities  mxx,  niyy,  mzz,  mxy,  myz,  mxz,  and 
equivalent  dipole  coordinates  x®,  yo,  and  zo,  are  given  by  the  diagonal  elements  of  the 
covariance  matrix 

cov(m)  =  (FtF)_1Ft  cov(d)  F  (FTF)''  .  (2.1.20) 

For  magnetic  field  measurements  with  squared  uncertainty  a  and  noise  uncorrelated 
between  receivers,  cov(d)  =  diag(a  )  is  a  diagonal  matrix,  and 

cov(m)  =  (FTF)_1a2  .  (2.1.21) 

For  data  with  unequal  uncertainties,  the  rows  of  equations  (2.1.10)  and  (2.1.18)  are 
normalized  by  dividing  by  the  corresponding  uncertainties,  giving  the  normalized  data 
unit  uncertainties. 
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2.1.5  Principal  Moment  and  Principal  Direction  Uncertainties 

The  leading  six  by  six  sub-matrix  of  cov(  m )  gives  the  covariance  of  the  non- 
redundant  elements  of  the  dipole  polarizablity  matrix  M,  cov(m).  The  principal 
directions  of  M  are  given  by  the  eigenvectors  of  M,  and  form  the  columns  of  the  rotation 
matrix  U  which  diagonalizes  M  (equation  2.1.7),  yielding  its  the  principal  moments  on 
the  diagonal.  Using  the  symmetry  of  M  and  L,  equation  (2.1.7)  can  be  written  as 

lL  =  Om  (2.1.22) 

where  1l  =  (Ln,  L22,  L33,  L12,  L23,  Lb)t,  and  O  is  an  orthogonal  matrix  obtained  by 
writing  out  matrix  product  (2.1.7)  explicitly  and  identifying  coefficients.  Principal 
moments  Ln,  L22,  L33  are  Rayleigh  quotients  of  matrix  M,  so  are  insensitive  to  first  order 
to  changes  in  estimated  principal  direction  matrix  U.  Their  squared  uncertainties  lie  on 
the  diagonal  of 

cov  (1L)  =  O  cov  (m)  Ot  .  (2.1.23) 

Uncertainties  in  the  principal  directions  of  M  are  related  to  the  stability  of  eigenvectors  of 
M  to  changes  in  M.  Perturbing  M  by  AM,  the  resulting  change  in  the  j'th  eigenvector 
(principal  direction)  Uj  is 


k*j 


uk  AM  Uj 


(2.1.24) 


to  first  order  in  AM,  provided  that  Xj  *  Xk  for  k  ^  j,  where  Xk  are  eigenvalues  of  M  (Ln, 
L22,  and  L33)  (Watson,  1983).  The  numerator  can  be  written  as 

ukT  AM  Uj=wjkTAm  (2.1.25) 


where  wjk  =  (uyuik,  u2jU2k,  u3ju3k,  uijU2k  +  u2jUik,  u2jU3k  +  u3jU2k,  uijU3k  +  u3juik).  The 
squared  uncertainties  of  the  elements  of  the  j'th  principal  direction  are  then  given  by  the 
diagonal  elements  of  the  three  by  three  matrix 

\T 


COV(Uj)  : 


(  T ) 

(  j  \ 

X-  UkWjk 

cov(m) 

^  «kWjk 

44 

1 

(2.1.26) 


If  some  eigenvalue  Xk  is  very  close  to  X}  the  denominator  in  equation  (2.1.24)  becomes 
small,  and  a  perturbation  AM  may  perturb  the  j'th  eigenvector  a  large  amount  in  the 
direction  of  the  k'th  eigenvector.  Consequently,  principal  directions  corresponding  to 
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two  principal  moments  become  indeterminate  when  the  difference  between  the  two 
moments  is  less  than  the  uncertainty  in  the  difference.  The  squared  uncertainty  in  the 
difference  between  the  i'th  and  j'th  principal  moments  is 

var  (Lii  -  Ljj)  =  cov  (lL)ii  +  cov  (lL)jj  -  2  cov  (lL)ij  ,  (2. 1 .27) 

where  cov  (lijij  is  the  i  j'th  element  of  cov  (1l). 

2.1.6  Application 

Our  current  application  of  equivalent  dipole  polarizabilities  is  discrimination 
amongst  buried  metallic  objects.  The  authors’  encoding  of  the  preceeding  algorithms 
have  been  extensively  tested  on  synthetic  data.  Two  synthetic  examples  are  presented 
here. 

The  first  example  simulates  collection  of  magnetic  induction  data  in  the  vicinity 
of  a  12  cm  diameter  buried  steel  sphere  with  a  relative  permeability  p,=l  80,  and 
conductivity  a  =107  Q"1  m'1,  with  the  sphere  center  1  m  below  the  level  of  transmitter  and 
receiver  coils.  Three  components  of  the  time  derivative  of  the  secondary  magnetic 
induction  dB(s)/dt  were  computed  at  the  center  of  a  1  m  square  loop  transmitter  for  81 
placements  of  the  loop  on  a  9  x  9  grid  with  0.4  m  spacing.  An  observation  time  of  610  ps 
after  transmitter  turn-off  was  chosen  to  approximate  the  effective  center  time  of  the 
averaging  gate  of  a  commercial  transmitter-receiver  system  (Geonics  EM-61).  The 
largest  observed  derivative  component  is  dBz(s)/dt  directly  above  the  sphere.  For  a  180 
Amp-m2  transmitter  moment,  dBzls)/dt  =  -4648.  nT/s  for  the  measurement  directly  above 
the  sphere  at  610  ps.  Gaussian  noise  of  magnitude  8.8  nT/s  was  added  to  the  dBz(s)/dt 
measurements  simulating  an  observed  noise  level  (at  Fort  Ord,  California).  Gaussian 
noise  of  magnitude  27.  nT/s  was  added  to  the  dBx(s)/dt  and  dBy(s)/dt  measurements  to 
simulate  the  larger  noise  levels  typically  observed  in  horizontal  field  components.  These 
data  were  inverted  for  dipole  polarizabilities  and  location.  The  downhill  simplex 
algorithm  converges  to  a  weighted  rms  misfit  of  0.90318  with  the  estimated  object  center 
at  (x,y,z)  =  (0.0028,  -0.0039,  1.0008)  meters.  Started  from  this  point,  after  three 
iterations  the  linearized  inversion  converges  to  a  weight  rms  misfit  of  0.90316  with  the 
estimated  object  center  at  (0.0026  ±  0.0030,  -0.0040  ±  0.0030,  1.0002  ±  0.0051)  meters. 
The  true  center  position  is  (0,0,1)  meters.  The  estimated  principal  dipole  polarizabilities 
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Ln,  L22,  and  L33  are  -0.654  ±  0.015,  -0.647  ±  0.011,  -0.635  ±  0.014  Amp-m2  /  s  /  pT. 
The  absolute  differences  between  principal  dipole  polarizability  estimates  0.008  ±0.010 
and  0.012  ±  0.009  Amp-m2  /  s  /  pT  respectively  for  |Ln  -  L22I  and  IL22  -  L33I,  are  less  than 
two  estimation  errors,  indicating  that  the  object  is  spherically  symmetric  within 
measurement  errors. 

For  a  second  example,  the  response  of  an  aluminum  prolate  spheroid  24  cm  long 
by  8  cm  wide,  of  conductivity  a  =  3.5  107  Q'1  m'1  was  modeled  using  an  integral  equation 
code  provided  by  P.  B.  Weichman  of  Blackhawk  Geophysics,  with  subsequent 
modifications  to  improve  accuracy.  The  code  expands  the  electric  field  within  the 
spheroid  in  a  polynomial  basis,  and  solves  for  a  set  of  modes,  each  with  a  characteristic 
decay  time.  Subsequently,  the  excitation  of  the  modes  for  each  position  of  transmitter 
loop  is  computed  for  a  ramp-on/ramp-off  transmitter  current,  and  the  contributions  of  the 
different  mode  voltages  observed  in  receiver  coils  are  summed  over  modes,  for  each 
transmitter-receiver  pair.  This  code  was  used  to  compute  the  spheroid  response  for  81 
transmitter  loop  positions  of  a  1  m  square  horizontal  loop,  on  a  9  by  9  grid  with  0.2  m 
spacing,  in  two  coaxial  dipole  receivers,  one  concentric  with  the  transmitter,  and  the  other 
0.4  m  above  the  first.  The  spheroid  center  was  placed  0.6  m  below  the  transmitter  level, 
offset  0.2  m  in  x  and  y  from  the  grid  center,  with  symmetry  axis  in  the  y-z  plane  dipping  - 

O 

30  .  A  3.3  ms  ramp-on,  0.08  ms  ramp-off  transmitter  current,  and  a  0.4  ms  averaging 
gate  starting  0.42  ms  after  transmitter  current  extinction,  were  used  to  emulate  a 
commercial  transmitter-receiver  system  (Geonics  EM-61). 

Gaussian  noise  with  a  magnitude  of  1%  of  the  largest  observed  voltage  was  added 
to  the  computed  voltages.  The  resultant  data  was  inverted,  yielding  an  estimated  center 
location  of  (0.207  ±  0.008,  0.206  ±  0.009,  0.600  ±  0.005)  meters,  in  agreement  with  the 
true  location  (0.200,  0.200,  0.600).  The  principal  polarizabilities  were  estimated  as  0.785 
±  0.023,  0.768  ±  0.018,  and  0.529  ±  0.011  V/pT,  with  differences  0.016  ±  0.025  and 
0.238  ±  0.020  V/pT.  The  agreement  of  Ln  and  L  22  indicates  an  object  that,  within 
measurement  errors,  is  rotationally  symmetric  about  the  third  principal  direction.  The 
differences  between  the  third  moment  and  the  other  two  are  well  resolved,  and,  for  a  non¬ 
magnetic  object,  consistent  with  the  smaller  cross  section  perpendicular  to  the  symmetry 
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axis.  The  third  principal  direction  is  estimated  as  (0.016  ±  0.023,  0.850  ±  0.026,  -0.526  ± 
0.043),  in  agreement  with  the  true  axis  of  symmetry  (0,  0.866,  -0.500). 

For  the  transmitter-receiver  configurations  used  in  the  above  examples  (which 
each  have  more  than  one  receiver  component  of  data),  examination  of  plots  of  the 
squared  data  misfit  y2  as  a  function  of  candidate  center  position  (not  shown),  on  sections 
through,  and  near,  the  true  object  location,  typically  show  a  large  "valley"  in  data  misfit 
sloping  downwards  to  a  well  defined  minimum  at  the  object  center  (within  measurement 
errors),  with  small  shallow  secondary  valleys  close  to  the  transmitters  and  receivers. 
Close  to  the  transmitters  and  receivers,  misfit  topography  is  on  the  scale  of  the  data  grid 
spacing.  The  existence  of  such  topography  recommends  measurement  grid  spacings  finer 
than  the  height  of  the  transmitters  and  receivers  above  the  shallowest  expected  depth  of 
objects  (e.g.,  ground  surface). 

For  the  commonly  occurring  case  of  a  transmitter-receiver  system  with  a  single 
horizontal  transmitter  loop  and  single  coaxial  receiver,  plots  of  squared  data  misfit  as  a 
function  of  candidate  center  position  are  more  problematic.  Figure  2.1.1  shows  a  detail 
of  such  a  plot,  for  a  1  m  horizontal  loop  transmitter,  concentric  Bz  receiver  system,  sited 
on  a  9  by  9  grid  centered  over  a  sphere  at  (0,0,1)  meters.  There  is  a  clear  minimum  at  the 
true  object  center  (0,0,1),  but  also  a  shallow  local  minimum  at  (0.00,  0.00,  1.08)  meters. 
The  existence  of  local  minima  close  to,  but  distinct  from,  the  minimum  associated  with 
the  true  object  position,  means  that  one  must  be  extremely  cautious  in  interpreting  data 
from  single  transmitter  single  receiver  systems.  For  comparison,  Figure  2.1.2  shows 
squared  data  misfit  y2  over  the  same  region  of  candidate  object  center  position,  for  a 
similar  system  with  an  additional  Bz  receiver  0.4  m  above  the  first.  The  added  data 
eliminates  the  secondary  local  minimum. 
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weighted  squared  misfit  as  a  function  of  trial  object  position 
Vertical  section,  y  =0.000 

Figure  2.1.1:  Squared  data  misfit  x2  as  a  function  of  candidate  object  center  position  ro,  for 
1  m2  square  transmitter  loop,  concentric  vertical  dipole  receiver  system  sited  on  a  9  x 
9  grid,  centered  1  m  above  a  12  cm  steel  sphere  (detail). 


weighted  squared  misfit  as  a  function  of  trial  object  posilion 
Vertical  section,  y  =  0.000 


Figure  2.1.2:  Squared  data  misfit  x2  as  a  function  of  candidate  object  center  position  ro,  for 
1  m2  square  transmitter  loop,  2  coaxial  vertical  dipole  receiver  system  sited  on  a  9  x  9 
grid,  centered  1  m  above  a  12  cm  steel  sphere  (detail). 
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2.2  Depths  of  equivalent  dipole  polarizability  resolution  for  some  transmitter  receiver 
configurations 

Equivalent  dipole  polarizability  matrices  and  dipole  locations  are  a  convenient 
way  to  summarize  active  source  induced  current  magnetic  field  measurements  in  an 
interpretable  form.  The  matrices'  principal  moments  give  information  on  the  rotational 
symmetry  of  a  conductive  object,  their  principal  directions  yield  the  object's  orientation, 
and  the  dipole  location  ro  estimates  the  object's  center  position.  An  algorithm  for 
estimating  these  and  their  uncertainties  was  outlined  in  a  previous  section.  Here 
uncertainty  estimates  are  used  to  compute  the  depths  to  which  the  polarizability  matrices 
and  dipole  locations  can  be  estimated  for  steel  spheres  of  varying  radius,  for  several 
transmitter-receiver  configurations . 

Equivalent  dipole  polarizability  matrices  M  model  observed  secondary  magnetic 
fields  B(s)  (r,t),  in  terms  of  the  magnetic  fields  of  unit  dipoles  in  the  x,  y,  and  z 
directions,  Bx(d)(r),  Byld,(r),  B7(d|(r)  centered  at  some  location  ro,  and  the  primary 
(inducing)  magnetic  field  strength  Bl0)-g(t),  at  ro,  for  a  given  time  variation  g(t)  of 
primary  magnetic  field: 

B(s)  (r,  t)  =  I" Bx(d)(r),  By(d)(r),  Bz(d)(r)l  M(t)  B(0)  (2.2.1) 

In  this  model,  the  polarizability  matrix  is  independent  of  transmitter  and  receiver  geometry 
and  object  location,  depending  only  on  the  innate  properties  of  the  object,  its  orientation, 
and  the  transmitter  waveform.  The  principal  values  of  M,  depend  only  on  the  object, 
independent  of  its  orientation,  and  on  the  transmitter  waveform.  Equivalent  dipole  position 
ro  is  generally  assumed  to  coincide  with  the  object  center. 

For  typical  time  domain  systems,  secondary  fields  are  measured  after  primary  fields 
are  extinct,  at  which  time  the  entire  magnetic  field  is  secondary.  In  simplest  form,  the  time 
variation  of  the  primary  field  is  incorporated  in  the  estimated  polarizability  matrix  M(t), 
and  the  primary  magnetic  field  strength  vector  at  object  center  B(0),  for  a  given  source,  is 
simply  the  magnetic  field  there  for  a  D.C.  current  in  the  transmitter  coils  of  the  transmitter's 
nominal  current  strength  (e.g.,  peak  current  strength). 
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The  methods  of  Section  2.1  were  used  to  estimate  equivalent  dipole  polarizabilities 
from  synthetic  three  component  magnetic  field  data  for  a  vertical  magnetic  dipole  source  at 
13  sites  placed  symmetrically  on  a  grid  with  0.4  m  spacings  in  x  and  y  centered  1  meter 
above  steel  spheres  of  varying  radius,  modeled  with  a  conductivity  of  a  =  107  Q'1  nf 1  and 
relative  permeability  p,  =  180.  A  step  function  turn-off  transmitter  current  was  used,  as  the 
most  generic  of  waveforms,  and  an  observation  time  of  610  ps  after  turnoff  chosen  to 
simulate  the  effective  center  time  of  the  averaging  gate  of  an  existent  commercial 
transmitter-receiver  system.  Polarizability  estimates  are  listed  in  Table  I.  In  principle,  for 
spherically  symmetric  objects,  the  three  principal  moments  are  identical.  The  variation 
between  estimated  principal  moments  is  less  than  1%  for  spheres  smaller  than  25  cm 
radius.  The  variations  are  due  to  the  limited  spatial  extent  of  the  data  used,  the  presence  of 
non-dipole  moment  components  in  the  data,  and  truncation  of  the  data  at  four  significant 
figures.  In  subsequent  computations  based  on  Table  I,  the  three  estimated  moments  were 
replaced  with  their  average. 


Radius  (m) 

Estimated  Principal  Moments  (Amp-m2/s/T) 

0.01 

-7.283  103 

-7.275  103 

-7.275  103 

0.02 

-1.114  104 

-1.114  104 

-1.113  104 

0.03 

-5.262  104 

-5.262  104 

-5.262  104 

0.04 

-1.530  105 

-1.530  105 

-1.529  105 

0.05 

-3.405  105 

-3.405  105 

-3.399  105 

0.06 

-6.416  105 

-6.416  105 

-6.401  105 

0.08 

-1.672  106 

-1.672  106 

-1.664  106 

0.10 

-3.382  106 

-3.382  106 

-3.355  106 

0.15 

-1.121  107 

-1.121  107 

-1.112  107 

0.25 

-4.418  107 

-4.418  107 

-4.373  107 

0.50 

-2.353  108 

-2.353  108 

-2.285  108 

1.00 

-1.098  109 

-1.098  109 

-1.057  109 

Table  I.  Estimated  principal  moments  for  simulated  steel  sphere  data  at  610  ps  after  step 
function  turn-off,  from  13  vector  measurements  coincident  with  vertical  dipole  source  on 
grid  1  m  above  sphere  center. 
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Section  2.1  gives  equations  for  the  covariance  matrix  for  dipole  polarizability 
matrix  M  elements,  and  covariance  matrices  for  the  principal  moments  and  directions 
derived  from  it.  In  general,  computing  principal  moments  and  directions  requires  knowing 
all  elements  of  M  (which  is  symmetric).  If  principal  directions  are  known  a  priori , 
principal  moments  may  be  determined  from  fewer  measurements,  but  determining  the 
principal  directions  empirically  requires  knowledge  of  all  elements  of  M.  A  simple 
measure  of  how  well  a  data  set  resolves  M  is  the  relative  average  squared  uncertainty 

^  =  Z  X  var  (mij )  1  Z  X  mij2  (2-2-2) 

i=l  j=l  i=l  j=l 

where  var(niij)  is  the  estimated  squared  uncertainty  of  the  ij'th  element  of  M  obtained  from 
the  diagonal  of  the  covariance  matrix  of  the  elements  of  M  (cov(m  )  of  Section  2.1).  The 
denominator  in  (2.2.2)  is  a  matrix  invariant;  independent  of  the  coordinates  used  to  express 
M,  so  a  property  of  the  object  from  which  M  arises.  In  our  experience,  the  numerator  of 
(2.2.2)  is  also  coordinate  independent,  suggesting  that  q2  itself  may  be  a  coordinate 
independent  measure  of  the  relative  uncertainty  in  M. 

For  each  radius  sphere  in  Table  I,  for  each  of  a  number  of  transmitter-receiver 
configurations,  the  relative  root  mean  squared  (rms)  moment  uncertainty  \  was  computed 
as  a  function  of  sphere  depth,  for  spheres  directly  below  the  center  of  a  9  x  9  grid  of  system 
placements  with  0.4  m  spacing  in  x  and  y.  One  meter  square  transmitter  loops  were  used 
with  a  moment  of  180  Amp-m2,  and  a  receiver  noise  level  of  1.97  nT/s  in  vertical  field 
measurements,  simulating  an  observed  noise  level,  and  5.91  nT/s  in  horizontal  field 
components  (when  present)  simulating  the  larger  noise  levels  observed  in  horizontal 
components. 

A  plot  of  relative  rms  moment  uncertainty  £,  as  a  function  of  depth  below 
transmitter  and  receiver  is  shown  in  Figure  2.2.1  for  an  isotropic  -6.41  10  Amp-m  /s/T 
equivalent  dipole  polarizability  (6  cm  radius  steel  sphere  at  610  ps)  below  the  grid  of 
system  placements  for  a  horizontal  loop  transmitter  /  concentric  vertical  magnetic  dipole 
receiver  system.  Being  based  on  a  linearized  inversion  for  M  and  ro,  these  uncertainty 
estimates  scale  linearly  with  receiver  noise  level.  For  spheres  very  near  the  level  of  the 
transmitter  and  receiver  (z=0),  the  rms  uncertainty  is  large  as  all  the  transmitter  placements 
illuminate  the  sphere  with  nearly  vertical  primary  fields,  yielding  little  information  on  mxx, 
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mxy,  and  myy,  and  correspondingly  large  variances  in  them.  Relative  rms  uncertainty 
decreases  to  a  minimum  at  0.135  m  depth,  and  subsequently  rises  with  increasing  depth, 
reflecting  the  decrease  of  primary  field  strengths  with  increasing  depth.  In  practice, 
mounting  a  receiver-transmitter  system  above  ground  level  sets  a  minimum  depth  below 
transmitter  for  buried  objects,  avoiding  difficulties  with  the  large  uncertainties  at  the 
transmitter-receiver  level.  The  great  increase  of  uncertainty  with  depth  for  spheres  at  large 
depths  limits  the  depths  for  which  polarizabilities  can  be  resolved.  For  the  horizontal  loop 
transmitter  vertical  dipole  receiver  on  this  grid  E,  >  0.1  for  6  cm  radius  steel  spheres  below 
1.47  m.  (At  z=1.47  m,  the  corresponding  relative  uncertainties  in  horizontal  and  vertical 
polarizabilities  are  0.061,  0.061,  and  0.147  respectively.) 


depth  (m) 

Figure  2.2.1:  Relative  rms  polarizability  uncertainty  E,  as  a  function  of  sphere  center  depth 
below  transmitter  and  receiver  for  a  6  cm  radius  steel  sphere.  Also,  relative 
uncertainty  in  vertical  moment  dmzz/dt  (dashed),  and  in  horizontal  moments  dmxx/dt 
and  dniyy/dt  (dotted). 

The  depths  to  5,  10,  and  20%  rms  uncertainty  in  polarizability  ds%(p),  dio%(p),  and  d2o%(p), 
were  similarly  found  for  isotropic  polarizabilities  corresponding  to  all  the  spheres  in  Table 
I,  and  are  plotted  in  Figure  2.2.2,  as  a  function  of  sphere  radius.  The  depths  to  5,  10,  and 
20%  uncertainty  in  estimated  sphere  depth,  &s%  z\  dio%  ,  and  d2o%  ,  are  plotted  in  Figure 
2.2.3.  In  general,  object  position  can  be  estimated  more  precisely  than  the  full 
polarizability  matrix  can,  as  object  position  may  be  determined  when  an  object  is 
illuminated  by  only  a  single  orientation  of  primary  field,  whereas  estimating  the  full 
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polarizability  matrix  requires  illuminating  the  object  with  primary  fields  Bo  in  at  least  three 
directions,  each  with  a  significant  component  in  the  direction  orthogonal  to  the  other  two. 
Consequently,  object  depth  can  be  resolved  within  10%  to  greater  depths  than  polarizability 
in  all  cases  plotted.  The  relative  uncertainties  in  polarizability  and  position  depend  on  all 
elements  of  the  polarizability  matrix  M  (through  F,  g,lq),  and  nik/q)  of  Section  2.1), 
depending  on  both  its  principal  values,  and  principal  directions  (object  orientation).  Most 
cases  presented  here  are  for  spherical  objects,  for  which  the  principal  values  are  all  the 
same,  and  the  polarizability  M  is  coordinate  independent.  For  comparison  dio%(p)  and 
dio%<z)  are  plotted  in  Figure  2.2.4  for  the  same  transmitter-receiver  pair,  for  the  case  of 
objects  with  the  same  vertical  dipole  polarizabilities  d  mzz  /dt  at  610  ps  as  the  spheres  of 
Table  I,  and  all  other  polarizabilities  null,  corresponding  to  thin  horizontal  non-magnetic 
discs.  The  general  trends  are  the  same  as  for  the  sphere.  Polarizability  can  be  resolved  to 
10%  slightly  deeper  than  for  the  sphere  for  all  but  the  1  cm  radius  sphere.  Object  depth  can 
be  resolved  approximately  1 .2  times  deeper  than  for  the  sphere. 


Depth  to  5%,  10%,  20%  moment  uncertainty 


radius  (m) 

Figure  2.2.2:  Depths  to  5%,  10%,  and  20%  rms  uncertainty  in  polarizability  as  a  function 
of  sphere  radius  for  steel  spheres  below  a  9  x  9  grid  of  transmitter-receiver  system 
placements,  for  aim2  square  loop  transmitter  with  a  concentric  vertical  dipole 
receiver. 
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Depth  to  5%,  10%,  20%  uncertainty  in  z 


radius  (m) 


Figure  2.2.3:  Depths  to  5%,  10%,  and  20%  uncertainty  in  sphere  center  depth,  as  a 

function  of  sphere  center  depth  below  center  of  transmitter-receiver  placement  grid, 
for  same  system  and  grid  as  in  Figure  2.2.2. 
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Figure  2.2.4:  Depths  to  10%  rms  polarizability  uncertainty  and  to  10%  uncertainty  in 
center  depth,  for  thin  horizontal  non-magnetic  disk  with  same  vertical  polarizability 
dm^/dt  at  610  ps,  as  spheres  of  Table  I,  plotted  as  a  function  of  the  corresponding 
sphere's  radius. 
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Effects  of  adding  a  second  coaxial  vertical  receiver  0.2  m  above  the  first  are  shown 
in  Figure  2.2.5,  where  dio%(p)  and  dio%<z)  are  plotted  for  the  two  receiver  system  (dotted) 
together  with  their  values  for  the  previously  plotted  single  receiver  system  (solid).  In 
addition  to  resolving  polarizabilities  and  location  to  greater  depth  as  shown  here,  the  added 
receiver  makes  locating  the  object  position  an  easier  problem  as  it  eliminates  a  secondary 
minimum  in  data  misfit  near  the  true  object  position  (Section  2.1). 


Depth  to  10%  uncertainty 


radius  (m) 

Figure  2.2.5:  Depths  to  10%  rms  polarizability  uncertainty  and  to  10%  uncertainty  in 
center  depth,  as  a  function  of  radius  for  steel  spheres  below  aim2  square  transmitter 
loop  with  two  coaxial  vertical  component  receivers  0.2  m  apart  vertically,  on  same 
grid  as  in  Figure  2.2.2. 

Figure  2.2.6  plots  where  dio%(p)  and  dio%(7)  (both  dotted)  for  a  similar  system  with 
the  two  vertical  component  receivers  in  the  plane  of  the  transmitter  loop,  offset  ±  0.2  m  in  x 
and  y  along  a  diagonal  from  the  loop  center.  This  system  shows  greatest  improvement  in 
sensitivity  over  the  single  receiver  system  for  objects  close  to  the  transmitter-receiver 
plane,  greatly  increasing  the  depths  to  which  the  smallest  spheres  can  be  resolved.  For 
spheres  of  radii  greater  than  8  cm,  the  (modest)  improvement  is  due  primarily  to  a  simple 
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reduction  in  uncertainty  by  a  factor  of  1/V2  due  to  the  doubling  of  the  number  of 
receivers. 


Depth  to  10%  uncertainty 


radius  (m) 


Figure  2.2.6:  Depths  to  10%  polarizability  uncertainty  and  to  10%  uncertainty  in  center 
depth,  as  a  function  of  sphere  radius,  for  1  m2  loop  transmitter  system  with  two 
vertical  component  receivers  0.566  m  apart  on  diagonal  in  plane  of  transmitter  loop, 
on  same  grid  as  in  Figure  2.2.2. 


Results  for  a  horizontal  loop  transmitter  three  component  concentric  receiver  are 
shown  in  Figure  2.2.7  (dotted  lines).  For  comparison,  a  system  with  a  1  m2  horizontal  loop 
and  two  orthogonal  1  m  vertical  loop  transmitters  with  lower  edges  at  the  level  of  the 
horizontal  loop  is  shown  in  Figure  2.2.8  (dotted  lines).  The  results  for  the  3  transmitter  1 
receiver  system  are  substantially  better  than  for  the  1  transmitter  3  receiver  system,  in  part, 
because  of  the  greater  noise  level  in  horizontal  component  receivers. 
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Depth  to  10%  uncertainty 


radius  (m) 


Figure  2.2.7:  Depths  to  10%  polarizability  uncertainty  and  to  10%  uncertainty  in  depth  as 
a  function  of  sphere  radius,  for  1  m2  loop  transmitter  system  with  3  component 
concentric  receiver,  on  same  grid  as  in  Figure  2.2.2. 


Depth  to  10%  uncertainty 


radius  (m) 

Figure  2.2.8:  Depths  to  10%  polarizability  uncertainty  and  to  10%  uncertainty  in  depth,  as 
a  function  of  sphere  radius,  for  three  orthogonal  1  m2  loop  transmitter  system  with 
vertical  component  receiver  at  horizontal  loop  center,  on  same  grid  as  in  Figure  2.2.2. 
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As  a  final  example,  adding  2  orthogonal  horizontal  component  receivers  to  the  3 
transmitter  system  to  make  a  3  transmitter  3  receiver  system  yields  results  shown  in  Figure 
2.2.9.  The  added  horizontal  components  substantially  increase  the  depth  of  resolution  of  the 
1  cm  radius  sphere,  but  little  affect  results  for  3  cm  radius  and  larger  spheres. 


Depth  to  10%  uncertainty 


radius  (m) 


Figure  2.2.9:  Depths  to  10%  polarizability  uncertainty  and  to  10%  uncertainty  in  depth,  as 
a  function  of  sphere  radius,  for  3  transmitter  loop  system  with  3  component  receiver  at 
horizontal  loop  center,  on  same  grid  as  in  Figure  2.2.2. 
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3.  PERFORMANCE  ASSESSMENT  AND  DIRECT  OPTIMIZATION 


3.1  Introduction 

The  broad  definition  of  optimum  in  this  report  is  an  active  EM  system  that  can 
extract,  from  the  measurements,  the  best  possible  estimates  of  the  location,  orientation, 
depth,  size,  shape,  and  metal  content  of  a  buried  metallic  object  in  the  presence  of  the 
interfering  response  of  non-UXO  metallic  objects.  Discrimination  can  be  achieved  partly 
through  selective  filtering  of  the  response  inherent  in  the  system  design  and  partly  through 
post  acquisition  processing.  These  parameters  are  never  determined  perfectly  because  of  a 
fundamental  non-uniqueness  in  the  solutions  (the  response  of  the  target)  to  the  governing 
equations.  Also,  practical  limitations  introduced  by  depth  of  the  target,  the  response  of  a 
general,  inhomogeneous,  ground  and  the  presence  of  response  from  other  non-UXO 
metallic  objects,  as  well  as  the  signal  to  noise  limitations  imposed  by  weight  and  power 
considerations,  introduce  constraints  on  what  is  achievable.  An  optimum  system  is 
bounded  by  these  theoretical  and  practical  considerations  and  at  the  end  represents  a 
compromise  that  detects,  discriminates  and  classifies  to  an  agreed  upon  criterion  or 
specification. 

There  have  been  several  theoretical  and  numerical  studies  of  target  response  for 
new  UXO  EM  systems,  none  of  which  have  addressed  the  issue  of  optimization. 
Laboratory,  field  and  theoretical  analyses  of  the  frequency  domain  GEM-3  (e.g.  Won  et  al., 
1997)  and  of  the  EM-61  and  its  variants  EM61HH,  EM63,  etc.  (e.g.  McNeill  and  Bosnar, 
1996;  Khadr  et  al.,  1998;  Ware,  2000)  for  a  range  of  targets  have  demonstrated  some  of  the 
detection/classification  properties  and  potential  of  these  systems.  The  studies  focused  on 
the  response  of  a  particular  system  to  various  targets  but  not  on  the  broader  issue  of  what 
system  would  perform  best  over  the  targets.  Some  of  these  systems  rely  on  the  frequency 
domain  (spectral)  or  time  domain  transient  signal  characteristics  to  classify  the  target  (e.g. 
Keiswetter  et  al.,  1999;  Barrow  et  al.  1996;  or  Ware,  2000).  Others  rely  on  both  the 
spectral  response  and  the  variation  of  response  with  transmitted  field  polarization  (e.g.  Bell 
et  al.,  1998;  Khadr  et  al.,  1998;  Pasion  and  Oldenburg,  2001;  Snyder  et  al.,  1999).  A 
fundamental  and  critical  technical  comparison  of  time  and  frequency  approaches  has  been 
tactfully  avoided. 
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Our  research  has  clearly  shown  the  importance  of  using  multiple  polarizations  of 
the  incident  magnetic  field  to  stimulate  the  principal  induced  magnetic  polarizations  of  the 
target.  The  second  step  in  the  identification/discrimination  process  is  then  to  measure  the 
resulting  secondary  fields  at  enough  points  in  space  to  uniquely  determine  these  principal 
polarizations  -  the  inversion  or  interpretation  process. 

[Fundamentally  the  incident  magnetic  field  induces  a  circulating  current  in  a  confined 
conductor.  This  current  produces  a  magnetic  dipole  moment,  usually  referred  to  as  M.  If 
this  moment  is  normalized  by  the  inducing  magnetic  field,  the  resulting  quantity  is  called 
the  magnetic  polarization.  We  use  both  terms  in  this  report.]  Finally  the  characterization 
depends  on  recovering  the  broadband  polarization  response. 

The  concept  of  identifying  an  object  through  the  measurement  of  secondary  fields 
arising  from  induced  magnetization  and  induced  currents  is  illustrated  in  the  following 
cartoon. 


Figure  3.1.1:  Induced  magnetization  and  currents  in  3D  bodies. 


In  any  conducting  permeable  body  placed  in  a  magnetic  field  there  are  two  types  of 
induced  magnetic  moments.  At  zero  frequency,  dc,  the  incident  magnetic  field  magnetizes 
the  body  so  that  it  acquires  a  static  magnetic  dipole  moment.  This  vector  moment  is  in  the 
direction  of  the  inducing  static  magnetic  field.  As  the  frequency  of  the  inducing  field  is 
increased  a  circulating  current  is  induced  to  flow  by  virtue  of  Faraday’s  Law.  This 
circulating  current  produces  a  magnetic  moment,  which  is  in  the  opposite  direction  to  the 
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inducing  magnetic  field.  The  net  magnetic  moment,  that  is  responsible  for  producing  the 
so  called  secondary  fields  measured  away  from  the  body,  is  a  function  of  the  frequency, 
conductivity,  permeability,  size  and  shape  of  the  body.  In  this  cartoon  the  induced 
principal  magnetic  moments  are  shown  in  black  and  the  induced  principal  electromagnetic 
moments,  for  a  particular  frequency,  caused  by  the  induced  currents  are  shown  in  red.  The 
characteristics  of  these  principal  moments,  and  how  they  change  with  frequency,  are  the 
basis  of  shape  and  metal  content  determination.  The  magnetic  and  EM  moments  are 
always  in  opposite  directions.  The  sphere  has  equal  moments  in  any  three  orthogonal 
directions.  The  disk  and  cylinder  have  two  identical  minor  axis  moments  and  a  different 
axial  moment,  and  the  flattened  ellipsoid  has  three  different  moments.  It  is  also  important 
to  see  from  this  cartoon  that  the  magnetic  moments  behave  differently  from  the  EM 
moments  and  can  be  large  in  directions  where  the  EM  moments  are  small.  For  example  in 
the  thin  disk  the  axial  magnetization  will  be  very  small  whereas  the  axial  moment  will  be 
the  largest  for  EM  induction.  [This  simplistic  description  for  the  disk  is  in  fact  only  valid 
for  very  thin  disks.  For  a  disk  of  appreciable  thickness  the  coplanar  EM  moments  can  be 
quite  large  due  to  the  increase  in  induction  number  caused  by  the  permeability.  A  relative 
permeability  of  200  effectively  increases  the  thickness  by  200  and  consequently  presents  a 
large  cross-sectional  area  to  the  inducing  electromagnetic  field.]  These  observations 
clearly  illustrate  the  need  for  broadband  coverage  not  only  to  separate  ferrous  from  non- 
ferrous  response  but  also  to  assist  in  relating  the  response  to  the  shape  of  the  body. 

For  regular  bodies  of  revolution  the  induced  moments  are  aligned  with  the 
symmetry  axes  of  the  object,  and  for  a  uniform  inducing  field  they  do  not  change  direction 
with  frequency.  For  an  irregular  object  such  as  twisted  scrap  metal,  the  moment  directions 
do  change  with  frequency.  For  actual  search  systems  the  dimensions  of  the  object  may  be 
comparable  to  the  distance  to,  and  size  of,  the  transmitter  and  the  moments  may  change 
size  and  direction  as  a  function  of  frequency.  In  any  event  it  is  these  induced  moment 
characteristics  that  allow  determination  of  the  object  parameters.  To  actually  excite  these 
moments  it  is  clear  that  several  polarizations  of  incident  field  are  required. 

A  rigorous  approach  to  the  interpretation  of  the  secondary  fields  measured  with  any 
particular  search  system  would  involve  an  inversion  scheme  to  determine  the  parameters  of 
an  object  that  best  reproduced  the  field  data  in  a  numerical  simulation  of  the  actual  system. 
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Forward  modeling  codes  for  solid  metallic  objects  of  arbitrary  shape  are  only  now 
becoming  available,  and  for  bodies  of  arbitrary  relative  permeability  and  arbitrary 
frequency  they  are  slow  even  on  very  fast  computers.  It  is  certainly  impractical  to  consider 
this  approach  for  real  time  processing  and  interpretation  in  the  field.  Further,  there  are  at 
the  moment  no  modeling  codes  for  ferrous  and  non-ferrous  shells  of  various  shapes.  Since 
these  are  needed  to  represent  typical  UXO  (there  being  little  interest  in  detecting  solid  iron 
cannon  balls)  another  approach  is  needed  to  represent  the  response  of  metallic  objects. 

For  purposes  of  design,  system  evaluation  and  interpretation  we  have  adopted  a 
simple  dipole  moment  representation  of  a  target.  This  is  basically  the  same  approach  used 
by  Khadr  et  al.  (1998),  Bell  et  al.  (1998),  Pasion  and  Oldenburg  (2001)  and  Baum  (1999). 
With  this  approximation  it  is  assumed  that  any  target  can  be  represented  by  three 
orthogonal  dipole  moments,  the  principal  dipole  moments  (PDM)  or,  if  normalized  by  the 
incident  field,  the  principal  dipole  polarizations  (PDP).  We  have  concluded  from  these 
studies  and  our  own  simulators  that  a  satisfactory  approach  to  characterizing  the  target  is  to 
recover  the  PDM’s  and  their  directions  as  a  function  of  frequency  or  time.  The  use  of 
PDM’s  greatly  simplifies  the  inversion  process  used  to  find  the  location  and  vector  PDM’s 
of  the  object.  Moreover  this  representation  is  a  practical  way  to  characterize  any  object 
since  it  basically  encapsulates  the  response  in  a  compact  manner.  A  rigorous  identification 
will  require  matching  of  the  PDM’s,  as  a  function  of  frequency  or  time,  to  a  catalogue  of 
objects  whose  PDM’s  have  been  predetermined. 

Smith  and  Morrison  (2004)  describe  an  inversion  algorithm,  which  locates  an 
object  and  recovers  its  PDM  for  a  specific  search  system  in  a  given  amount  of  system 
noise.  The  inversion  process  inherently  assigns  uncertainties  to  the  extracted  target 
parameters  and  location  based  on  the  uncertainties  (noise)  in  the  data.  These  uncertainties 
are  the  fundamental  data  used  to  assess  the  quality  of  the  response  and  to  estimate  ROC 
curves.  Smith  et  al.  (2004b)  used  this  inversion  algorithm  to  determine  the  depths  of 
detection,  to  a  given  uncertainty,  and  PDM’s  to  a  given  uncertainty,  of  simple  target 
spheres  for  several  search  configurations.  This  inversion  analysis  underlies  the  design 
methodology  described  in  this  section  of  the  report. 
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3.2  Optimizing  the  transmitter-receiver  configuration 


The  inversion  algorithm  described  in  Smith  and  Morrison  (2004),  and  used  for 
system  comparison  by  Smith  et  al.  (2004b)  is  the  fundamental  tool  used  in  this  more 
generalized  study.  In  a  sense  the  approach  used  here  can  also  be  cast  as  an  inverse 
problem:  to  determine  the  parameters  of  the  AEM  system  that  maximize  the  response  for  a 
given  target  subject  to  constraints  on  the  system  (size,  platform,  weight  and  power)  and  the 
anticipated  ambient  noise. 

As  with  a  rigorous  inversion,  the  parameters  of  a  proposed  system  are  varied 
incrementally  until  a  maximum  response  is  obtained  or  until  a  maximum  in  the  signal  to 
noise  ratio  is  obtained.  For  a  given  model  and  background,  and  with  given  system  and 
ambient  noise,  the  code  is  run  repeatedly  for  changes  in  the  system  parameters  until  the 
variances  in  the  estimated  target  parameters  are  minimized.  This  process  provides  an 
objective  method  for  finding  the  optimum  array  configuration  (i.e.  the  number  of 
components  and  their  spatial  disposition,  the  transmitter  moments  and  the  bandwidth) 
needed  to  obtain  the  best  estimate  of  the  depth  and  vector  principal  moments  of  the  target. 

The  first  step  in  such  a  process  can  be  seen  in  one  of  the  examples  shown  in  Smith 
et  al.  (2004b)  (Figure  3.2.1).  Here  uncertainty  estimates  are  used  to  compute  the  depths  to 
which  the  PDM’s  and  dipole  locations  can  be  estimated  for  steel  spheres  of  varying  radius, 
for  two  transmitter-receiver  configurations.  One  meter  square  transmitter  loops  were  used 
with  a  moment  of  180  Amp-m  ,  and  a  receiver  noise  level  of  1.97  nT/s  in  vertical  field 
measurements,  simulating  an  observed  noise  level,  and  5.91  nT/s  in  horizontal  field 
components  (when  present)  simulating  the  larger  noise  levels  observed  in  horizontal 
components.  A  step  function  turn-off  transmitter  current  was  used,  as  the  most  generic  of 
waveforms,  and  an  observation  time  of  610  ps  after  turnoff  chosen  to  simulate  the  effective 
center  time  of  the  averaging  gate  of  an  existent  commercial  transmitter-receiver  (T-R) 
system.  For  each  radius  sphere  and  for  each  T-R  configuration  the  relative  root  mean 
squared  (rms)  inverted  moment  uncertainty  and  depth  uncertainty  were  computed  as  a 
function  of  sphere  depth,  for  spheres  directly  below  the  center  of  a  9  x  9  grid  of  system 
placements  with  0.4  m  spacing  in  x  and  y. 
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Depth  to  10%  uncertainty 


Depth  to  10%  uncertainty 


Figure  3.2.1:  Depth  to  10%  polarizability  uncertainty  and  10%  uncertainty  in  depth  as  a 
function  of  sphere  radius  for  TxTyTz-Bz  and  TxTyTz  -  BxByBz  systems. 


The  results,  in  Figure  3.2.1,  are  shown  for  an  orthogonal  three  loop  transmitter  and  vertical 
field,  Bz,  receiver  and  for  the  same  transmitter  with  three  orthogonal  receivers.  The  results 
for  the  standard  horizontal  loop  with  a  single  vertical  receiver  are  shown  for  reference  as 
the  solid  line.  The  important  role  of  multiple  field  polarizations  at  the  target  is  easily  seen 
in  this  figure.  The  depth  of  detection  for  a  10  cm  radius  sphere  almost  doubles  with  the 
three-component  transmitter.  But  relatively  little  is  gained  by  adding  a  triaxial  receiver. 

In  general  it  was  found  that  the  object  position  can  be  estimated  more  precisely  than 
the  PDM’s.  The  object  position  may  be  determined  with  only  a  single  orientation  of 
primary  field,  whereas  estimating  the  full  polarizability  matrix  requires  illuminating  the 
object  with  primary  fields  in  at  least  three  directions,  each  with  a  significant  component  in 
the  direction  orthogonal  to  the  other  two.  Consequently,  object  depth  can  be  resolved 
within  10%  to  greater  depths  than  PDM’s. 

The  recovery  of  the  PDM’s  requires  measurements  using  several  polarizations  of 
the  incident  field.  How  these  are  provided  is  a  function  of  the  spatial  deployment  of  a 
system  as  well  as  the  number  of  polarizations  provided  by  different  transmitters  in  the 
system.  If  a  single  horizontal  loop  source  is  moved  in  discrete  intervals  over  an  object  then 
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it  is  illuminated  with  different  polarizations  by  virtue  of  its  changing  position  in  the  source 
dipole  field.  This  is  illustrated  in  the  cartoon  of  Figure  3.2.2. 


Figure  3.2.2 

In  the  numerical  experiments  by  Smith  et  al.  (2004b)  referred  to  above,  the  ‘data’  were  in 
fact  assumed  to  have  been  taken  on  a  grid.  An  illustration  of  the  effectiveness  of  an 
elementary  system  employing  a  single  horizontal  loop  receiver  within  a  horizontal  loop 
transmitter  is  shown  in  Figure  3.2.3. 
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Figure  3.2.3:  Principal  polarizations  (PDM)  and  location  for  a  dipping  ellipsoid  on  9  x  9 
grid  using  a  single  horizontal  loop  receiver  within  a  horizontal  loop  transmitter. 
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The  target  is  a  dipping  ellipsoid  and  the  response  is  obtained  on  a  9  by  9  grid  above  it.  The 
table  shows  that  the  principal  polarizations  (called  PDM’s  in  this  figure)  and  location  are 
very  well  recovered.  In  this  inversion  instrumentation  noise  was  included  but  location 
errors  were  not.  The  accuracy  of  the  inversion  now  depends  on  accuracy  of  the  positioning 
on  the  grid.  We  have  addressed  the  role  of  positioning  error  in  this  project  and  some 
examples  are  presented  below.  For  rapid  field  surveys  absolute  positioning  to  the 
centimeter  level  might  be  expensive  and  difficult  but  relative  positioning  on  the  scale  of  the 
target  depth  might  be  easier  to  achieve.  This  suggests  designing  systems  which  are  as  close 
as  possible  to  stand  alone:  multiple  transmitters  and  a  number  of  spaced  apart  receivers  that 
locate  and  characterize  the  target  from  a  single  system  set-up  in  the  vicinity  of  the  target. 
Such  a  system  also  suggests  a  search  procedure  involving  two  modes  of  operation.  First, 
the  object  is  located  with  a  simple  configuration  and  low  power,  narrow  band,  mode. 
Once  located,  the  system  switches  to  broadband,  multicomponent  mode  and  the  accurate 
depths  and  PDM’s,  and  identification  through  the  look-up  catalog  of  object  parameters. 
Finally,  anticipating  the  results  presented  below,  a  system  comprising  multiple  receivers 
will  require  sensors  that  are  smaller  than  the  loop  receivers  currently  employed.  In  this 
report  we  analyze  a  variety  of  system  configurations  that  are  optimum  for  detecting  objects 
and  for  determining  their  PDM’s  while  minimizing  their  complexity  and  minimizing  the 
number  of  sensors  required. 

The  general  simulation/inversion  code  was  used  to  investigate  the  number  of 
receivers,  and  their  orientations,  required  for  an  optimum  measurement  of  the  principal 
polarizabilities  and  depths.  In  general,  it  was  found  that  the  object  position  can  be 
estimated  more  precisely  than  the  PDM’s  for  the  same  object  depth.  The  object  position 
may  be  determined  with  only  a  single  orientation  of  primary  field,  whereas  estimating  the 
full  polarizability  matrix  requires  illuminating  the  object  with  primary  fields  in  at  least 
three  directions,  each  with  a  significant  component  in  the  direction  orthogonal  to  the  other 
two.  Consequently,  object  depth  can  be  resolved  within  10%  to  greater  depths  than  PDM’s 
with  the  same  uncertainty.  The  following  brief  description  and  examples  are  taken  from  a 
comprehensive  analysis  of  this  problem  that  is  discussed  in  more  detail  in  the  paper  by 
Smith  et  al.  (2005)  that  has  been  submitted  for  publication. 
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In  this  analysis,  transmitter  systems  are  comprised  of  one  or  several  rectangular 
loops  of  fixed  size,  and  a  number  of  receiver  coils  approximated  as  point  measurements. 
Three  families  of  designs  are  presented:  a)  systems  for  use  on  a  2-D  grid  of  positions  with 
negligible  error  in  relative  instmment  location,  b)  systems  for  use  on  a  line  of  positions 
with  negligible  error  in  relative  instmment  location,  and  c)  systems  for  stand  alone  use, 
insensitive  to  instmment  positioning  errors. 

The  general  approach  is  summarized  with  reference  to  Figure  3.2.4  for  a  2D  grid  of 
data.  The  color  plot  shows  the  uncertainty  in  the  determination  of  the  principal 
polarizabilities  for  a  test  object  (in  this  case  a  sphere  for  ease  of  illustration)  of  various 
diameters.  The  data  for  the  inversion  are  the  responses  for  the  simple  in-loop  configuration 
(similar  to  the  commercial  EM  61  system)  at  81  positions  of  the  9  x  9  grid  shown  on  the 
right.  The  contours  are  drawn  at  locations  in  the  vertical  section  beneath  the  array  where 
the  uncertainty  in  polarizability  is  equal  to  the  polarizability  itself  -  in  effect  a  signal  to 
noise  ratio  of  one.  As  mentioned  above  the  depth  can  be  determined  with  any  specified 
uncertainty  at  greater  depth  than  the  polarizability.  Using  Figure  3.2.1  as  a  guide,  the  depth 
of  a  60  mm  sphere  can  be  determined  to  within  10%  at  a  depth  50%  greater  than  depth  at 
which  the  polarizability  is  determined  to  within  10%.  Consequently  the  contours  in  plots 
like  those  of  Figure  3.2.1  can  be  used  as  rough  measures  of  the  depth  of  detection  of  the 
relevant  object.  For  example  the  polarizability  of  a  40  mm  sphere  could  be  determined 
with  an  uncertainty  of  less  than  its  polarizability  anywhere  in  the  section  above  the  contour 
labeled  40  mm  but  the  depth  would  be  well  determined  below  this  contour.  The  40  mm 
sphere  would  consequently  be  well  located  anywhere  to  a  depth  of  about  one  meter  and  in  a 
swath  of  plus  or  minus  1.0  m  around  its  horizontal  position.  In  comparing  configurations 
for  a  specific  target,  the  ‘best’  are  those  for  which  the  contours  for  that  target’s  diameter  are 
deepest  in  the  section  and  flattest.  We  have  found  this  graphic  presentation  to  be  an 
excellent  way  to  illustrate  the  range  and  accuracy  of  a  particular  configuration  of 
transmitter  and  receiver. 
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Figure  3.2.4:  Rms  uncertainty  in  polarizability  as  a  function  of  position  on  9  x  9  grid 
using  simple  in-loop  configuration. 

The  improvement  in  detectability  that  accrues  from  using  multiple  transmitters  and 
receivers  is  evident  in  the  results  shown  in  Figure  3.2.5  for  a  configuration  with  3 
orthogonal  transmitters  and  3  orthogonal  receivers  deployed  on  the  same  9x9  grid.  The 
detectability  of  the  40  mm  sphere  is  pushed  down  to  1.5  m  and  the  80  mm  sphere, 
previously  detectable  to  1.5m  is  now  pushed  down,  off  the  graph,  to  more  than  2.0  m. 

Next,  in  Figure  3.2.6,  the  uncertainty  in  polarizability  is  plotted  for  a  simple  two 
transmitter  three  receiver  (vertical)  configuration  deployed  at  21  positions  along  a  line  over 
the  target.  In  this  and  subsequent  multi  receiver  configurations  the  inversion  code  was  used 
to  determine  the  location  of  the  receivers  that  optimized  the  resulting  uncertainty  plot.  We 
found  that  the  orientation  of  the  sensors  was  less  important  than  the  number  used  and  their 
spacing.  Since  we  have  observed  that  ambient  EM  noise  fields  are  predominantly  in  the 
horizontal  plane,  we  have  elected  to  use  vertical  receivers  in  the  array  studies.  Further,  to 
minimize  the  footprint  of  the  system  we  imposed  penalties  on  the  inversion  for  receivers 
outside  the  perimeter  of  the  horizontal  transmitter. 

The  results  for  the  line  system  indicate  that  the  depth  range  is  reduced  a  little  and 
the  swath  width  is  considerably  reduced  over  the  grid  system  of  Figure  3.2.5.  Nevertheless 
it  is  important  to  note  that  comparable  results  would  be  obtained  if  successive  lines  were 
run  with  line  spacing  of  0.5  m.  The  grid  after  all  could  be  considered  as  multiple  parallel 
lines  with  spacing  of  only  0.4  m.  If  survey  cost  is  proportional  to  the  number  of  readings 
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then  there  are  appreciable  savings  in  the  multi-element  line  profiling.  We  will  return  later 
to  the  role  of  position  uncertainty  in  this  assessment. 


Figure  3.2.5:  Rms  uncertainty  in  polarizability  as  a  function  of  position  on  9  x  9  grid 
using  3  orthogonal  transmitters  and  3  orthogonal  receivers. 


Figure  3.2.6:  Rms  uncertainty  in  polarizability  as  a  function  of  position  using  two 
transmitters  and  three  receivers  in  a  profile  mode. 


Finally  we  experimented  with  single  or  stand  alone  configurations  which  could 
determine  depth  and  polarizations  from  a  single  position  in  the  vicinity  of  the  target.  In  all 
these  systems  we  used  a  three  component  transmitter.  The  results  for  a  four  receiver 
configuration  are  shown  in  Figure  3.2.7.  The  depth  of  detection  is  reduced  somewhat  and 
the  pattern  is  narrowed  and  skewed.  Another  important  measure  of  the  design  process  is 
illustrated  for  this  array  in  the  plot  of  Figure  3.2.8.  Here  the  uncertainty  in  polarizability  is 
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plotted  as  a  function  of  orientation  of  an  elongate  target  at  a  depth  of  one  meter.  The  plot 
clearly  shows  that  there  are  ‘blind  spots’  for  this  configuration  -  orientations  for  which  the 
number  of  T-R  pairs  is  inadequate  to  determine  the  orientation.  Smith  et  al.  (2005)  have 
determined  theoretically  that  13  T-R  pairs  are  required  to  determine  the  polarizabilities  so 
the  12  pairs  are  in  fact  not  enough.  Adding  a  fifth  receiver,  in  the  optimized  five  element 
array  of  Figure  3.2.9,  makes  a  dramatic  improvement  in  the  orientation  determination. 
Even  at  a  depth  of  1.6  m  there  are  no  blind  spots  for  this  array  as  seen  in  Figure  3.2.10, 

although  there  is  only  a  modest  improvement  in  the  depth  of  detectability,  Figure  3.2.9. 
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Figure  3.2.7:  Rms  uncertainty  in  polarizability  as  a  function  of  position  using  three 
transmitters  and  four  receivers  in  a  stand  alone  mode. 
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Figure  3.2.8:  Rms  uncertainty  in  polarizability  as  a  function  of  orientation  of  an  elongate 
target  at  a  depth  of  1  m. 
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Figure  3.2.9:  Rms  uncertainty  in  polarizability  as  a  function  of  position  using  three 
transmitters  and  five  receivers  in  a  stand  alone  mode. 


Figure  3.2.10:  Rms  uncertainty  in  polarizability  as  a  function  of  orientation  of  a  target  at 
1.6  m  depth  using  three  transmitters  and  five  receivers  in  a  stand  alone  mode. 


Allowing  the  inversion  code  to  ‘choose’  the  location  of  the  vertical  receivers  results 
in  a  rather  irregular  pattern  as  shown  on  the  right  of  Figure  3.2.9  as  an  example.  From  a 
fabrication  viewpoint  it  might  be  simpler  to  arrange  the  sensors  in  a  regular  pattern  dictated 
by  the  frame  or  structure  of  the  transmitters.  We  tried  the  regular  pattern  of  Figure  3.2.1 1 
and  found  very  little  reduction  in  detectability  compared  to  the  optimum  pattern.  The 
conclusion  of  studies  such  as  this  is  that  optimization  is  not  a  strong  function  of  variation  of 
sensor  position  around  the  optimum  one. 
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Figure  3.2.11:  Rms  uncertainty  in  polarizability  as  a  function  of  position  using  three 
transmitters  and  five  receivers  in  a  regular  pattern  dictated  by  transmitter  frame. 

By  using  an  ingenious  combination  of  horizontal  loops  one  can  provide  both  a 
vertical  and  horizontal  polarization  of  the  incident  field  at  a  target  beneath  the  array.  In 
Figure  3.2.12  the  horizontal  loops  consist  of  four  independent  rectangular  loops  and  a 
redundancy  of  9  vertical  receivers.  When  two  adjacent  loops  are  energized  with  the  same, 
say  clockwise,  current  flow  they  sum  to  produce  an  equivalent  large  loop  of  the  same 
polarity  which  produces  a  primary  field  directed  vertically  beneath  the  center  of  the  loop. 
When  energized  with  opposing  current  flows  a  multipole  field  is  produced  with  field  lines 
that  are  horizontal  beneath  the  loops.  Similar  energization  of  the  orthogonal  pair  of  loops 
produces  the  requisite  third,  horizontal,  polarization  of  the  field  incident  on  the  target. 

The  results  (Figure  3.2.12),  even  with  9  receivers,  are  not  as  good  as  those  obtained 
with  the  transmitter  consisting  of  three  orthogonal  loops  and  5  receivers,  Figure  3.2.1 1. 

Finally,  to  summarize  these  illustrative  examples,  it  is  important  to  note  that  the 
depth  of  detectability  for  the  3  transmitter  5  receiver  stand  alone  system,  Figure  2.2.11,  is 
almost  as  good  as  the  3  transmitter  3  receiver  system  when  the  latter  is  deployed  at  21 
positions  along  a  line  or  profile.  It  is  also  almost  as  good  as  the  current  industry  standard 
EM61  system  when  the  latter  is  deployed  on  a  grid  of  81  points,  Figure  3.2.4.  On  the  other 
hand  it  is  not  nearly  as  good  as  the  3  transmitter  3  receiver  configuration  deployed  on  the 
same  grid,  Figure  3.2.5. 
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Figure  3.2.12:  Rms  uncertainty  in  polarizability  as  a  function  of  position  using  four 
independent  horizontal  loops  and  9  vertical  receivers. 

The  results  to  date  were  obtained  assuming  that  there  were  no  errors  in  the  position 
of  the  system  on  the  grid  or  along  the  line.  In  a  practical  survey  there  are  such  errors, 
largest  in  the  absolute  position  on  the  grid,  less  in  relative  position  along  a  line  segment  and 
zero  for  the  stand  alone  system  (the  latter  finds  the  polarizability  and  the  target  location 
relative  to  the  T-R  system  so  the  uncertainties  are  proportional  only  to  the  system  noise). 
Figures  3.2.13  and  3.2.14  display  the  relative  uncertainties  in  polarizability  and  depth 
respectively  as  a  function  of  the  position  error  for  the  three  deployment  arrays  shown  on 
the  right  of  the  figures.  The  statistics  were  created  by  repeated  inversions  with  random 
misplacements  of  each  measurement  point  by  the  indicated  position  error.  The  target  in 
this  experiment  was  a  horizontal  22  mm  shell  with  an  aspect  ratio  of  8:1  at  a  depth  of  0.75 
m.  The  relative  uncertainties  in  polarizability  and  depth  for  both  the  line  and  grid 
configurations  remain  less  than  that  of  the  stand  alone  system  (about  5%)  as  long  as  the 
position  uncertainty  is  less  than  about  2  cm.  Relative  GPS  positioning  is  at  best  a  few  cm 
so  without  expensive  and  logistically  difficult  laser  or  microwave  positioning  it  appears 
that  the  stand  alone  system  is  superior  to  line  and  grid  deployments.  But  the  stand  alone 
system  can  also  be  operated  in  line  or  grid  mode  the  results  of  which  would  be  significantly 
better  than  any  of  the  individual  schemes  discussed  till  now.  Since  the  stand  alone  system 
has  the  same  footprint  as  the  best  current  commercial  system,  the  field  operational  issues 
are  identical  to  the  current  systems  but  with  vastly  improved  detection  and  identification 
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properties.  Further  the  system  retains  excellent  detection  and  polarization  determination 
from  irregularly  located  single  positions,  which  may  be  all  that  is  available  in  terrain  with 
obstacles  to  regular  grids  or  profiles. 


Horizontal  22mm  8:1  aspect  ratio  ‘shell’  at  0.75  m  depth 
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Figure  3.2.13:  Relative  rms  polarizability  uncertainty  as  a  function  of  instrument  location 
error  for  22  mm  8: 1  aspect  ratio  shell  at  0.75  m  depth  at  610  ps  after  transmitter  shut¬ 
off. 


Horizontal  22mm  8:1  aspect  ratio  ‘shell'  at  0.75  m  depth 
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Figure  3.2.14:  Relative  rms  depth  uncertainty  as  a  function  of  instrument  location  error  for 
22  mm  8:1  aspect  ratio  shell  at  0.75  m  depth  at  610  ps  after  transmitter  shut-off. 
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Given  the  position  and  polarizability  matrix  of  an  object,  it  is  straightforward  to 
calculate  estimation  uncertainties  in  polarizabilities  and  coordinates  that  would  arise  from 
detection  with  given  transmitter  and  receiver  configurations  with  specified  source  current 
and  receiver  noise  levels.  One  can  compare  performance  of  prospective  systems  using  a 
simple  statistic  such  as  the  total  summed  squared  uncertainty  in  polarizability  estimates 
for  on  object  at  a  specific  position  relative  to  the  detection  system,  or  sum  the  squared 
polarizability  uncertainties  for  test  objects  at  a  range  of  positions  (‘control  points’)  to 
construct  a  somewhat  more  general  measure  of  performance.  Since  uncertainties 
typically  vary  greatly  with  the  depth  of  an  object  below  a  system,  it  makes  sense  to 
include  objects  at  number  of  depths  in  the  sum,  and  to  weight  them  differently  according 
to  depth,  for  example,  using  the  weight  function  w2(z)  =  (z/zmax) 7  which  results  in 
roughly  equal  contributions  from  shallow  and  deep  objects. 

As  conceiving  and  comparing  a  large  number  of  prospective  configurations  is  a 
potentially  burdensome  task,  we  have  automated  the  procedure,  in  effect,  using  general 
non-linear  optimization  methods  to  minimize  the  weighted  summed  squares  of 
polarizability  uncertainties  for  prospective  objects  over  a  range  of  depths  below  the 
system.  For  simplicity,  we  keep  the  transmitter  configuration  fixed;  for  example,  three 
orthogonal  1  m2  rectangular  loop  transmitters  abutting  the  horizontal  plane  (z  =  0),  with 
centers  at  horizontal  position  x  =  y  =  0,  and  allow  receiver  positions  (and  possibly 
orientations)  to  vary  for  a  fixed  number  of  receivers.  To  model  the  typically  larger 
amount  of  noise  in  horizontal  field  measurements  compared  to  vertical  field 
measurements  we  treat  receiver  noise  as  three  times  greater  in  horizontal  than  vertical 
receiver  orientations  with  a  sinusoidal  variation  in  between. 

In  general,  we  find  that  for  spherical  objects  (three  equal  principal  polarizabilities) 
receiver  position  and  orientation  optimization  results  in  all  receivers  being  oriented 
vertically.  Optimization  is  simplest  for  spherical  test  objects  as  spheres  are  rotationally 
invariant.  Polarizability  estimation  uncertainties  do  depend  on  the  polarizability  matrix. 
Computing  summed  squared  polarizability  uncertainty  for  elongate  test  objects  (with  only 
a  single  non-zero  principal  polarizability  in  their  limiting  case)  one  can  compare  for 
uncertainties  for  elongate  objects  of  any  orientation.  Doing  this  for  systems  with  four, 
five,  or  six  receivers  optimized  for  spherical  targets,  one  finds  that  four  receiver  systems 
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have  a  continuous  series  of  elongate  object  orientations  for  which  the  summed  squared 
uncertainty  is  extremely  large  (or  unbounded).  For  five  receiver  systems,  a  few  elongate 
object  orientations  remain  with  extremely  large  uncertainties.  And  for  six  receivers  the 
problem  is  greatly  mitigated,  but  still  one  finds,  for  example,  some  horizontal  orientations 
result  in  rms  polarizability  uncertainties  50  times  greater  than  for  the  same  elongate 
object  in  a  near  vertical  orientation. 

In  performing  the  same  sort  of  optimization  for  non-spherical  objects,  one  has  the 
difficulty  that  polarizability  uncertainty  depends  on  object  orientation.  A  quick  method 
to  partially  stem  this  difficulty  is  to  simply  sum  squared  polarizability  uncertainties  over 
the  test  object  in  a  number  of  standard  orientations  at  each  of  the  control  points  at  which 
uncertainties  are  computed  for  the  test  object,  for  example,  orienting  the  object  major 
axis  along  directions  corresponding  to  the  face  centers  and  comers  of  a  cube  centered  at 
the  origin.  For  example,  doing  this  in  a  six-receiver  optimization  for  elongate  objects 
over  a  range  of  test  depths  from  0.2  m  to  1.6  m,  reduces  the  variation  in  rms  polarizability 
uncertainty  with  orientation  for  objects  at  1  m  depth  from  the  factor  of  50  to  a  factor  of  5, 
with  a  6.5  reduction  of  the  largest  rms  polarizability  uncertainty.  This  analysis  is 
reported  in  detail  in  Smith  et  al.  (2005). 

Worst-case  optimization 

Summing  squared  uncertainties  over  estimates  for  a  number  of  orientations  of  an 
elongate  object  is  a  simple  ad  hoc  method  of  adapting  the  minimization  objective  to  find 
receiver  arrays  which  perform  better  for  elongate  objects  than  systems  which  have  only 
been  optimized  for  resolving  spheres.  Being  ad  hoc,  it  leaves  some  uncertainty  as  to 
whether  there  might  be  other  objects  with  polarizabilities  that  would  be  unresolvable, 
despite  having  sizable  polarizabilities. 

To  first  order  in  receiver  noise  level,  for  a  given  transmitter-receiver  system  and 
object  position  and  orientation,  the  polarizability  uncertainties  are  independent  of  the 
scale  of  the  polarizability  matrix  (Smith  and  Morrison,  2004);  the  polarizability 
uncertainties  for  an  object  with  small  polarizabilities  and  for  a  similar  object  with  large 
polarizabilities  with  the  same  orientation  and  the  same  ratios  of  principal  polarizabilities 
are  the  same.  Of  course,  for  an  object  with  small  polarizabilities,  the  relative 
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uncertainties  in  the  polarizabilities  are  much  greater  than  the  relative  uncertainties  for  the 
similar  object  with  larger  polarizabilities.  So,  in  optimizing  a  system,  one  need  not  worry 
about  the  scale  of  the  polarizability  matrix  for  which  it  is  optimized. 

Ostensibly,  one  might  like  to  find  systems  that  are  optimal  in  a  mini-max  sense: 
giving  the  minimum  summed  squared  polarizability  uncertainty  for  whatever 
polarizability  object  has  the  greatest  summed  squared  uncertainty.  In  general,  as  we 
consider  both  object  polarizability  and  position  as  parameters  to  be  estimated, 
uncertainties  in  one  and  the  other  are  not  independent,  so  a  mini-max  problem  for 
uncertainties  in  one  separate  from  the  other  is  not  convenient.  For  mini-max 
minimization,  instead  of  considering  the  summed  squared  uncertainty  in  all  components, 
it  is  much  more  convenient  to  worry  about  only  the  squared  uncertainty  in  the  component 
that  is  least  well  constrained  by  a  set  of  data.  In  summing  over  squared  uncertainties,  this 
component  would  make  the  greatest  contribution  to  the  summed  squared  uncertainty. 
Given  the  structure  of  the  sensitivity  matrix  (Jacobian)  of  data  with  respect  to 
polarizabilities  and  position,  as  it  depends  on  object  polarizability,  we  find  the  worst-case 
polarizability  matrix,  which  has  the  component  with  the  largest  uncertainty  of  all 
components  of  all  polarizability  matrices.  In  most  cases  examined,  the  worst  case 
polarizabilities  have  one  principal  polarizability  that  is  less  than  one  hundredth  of  the 
largest  principal  polarizability,  and  a  second  principal  polarizability  that  is  less  one  tenth 
of  the  largest  principal  polarizability. 

For  a  given  receiver-transmitter  system  and  object  location,  once  the  worst  case 
polarizability  is  in  hand,  it  is  straightforward  to  compute  its  summed  squared 
polarizability  uncertainty,  and  to  use  this  summed  over  prospective  object  positions 
(control  points)  as  an  objective  function  in  non-linear  optimization.  More  details  are 
given  in  Smith  and  Morrison  (2005b). 

3.3  Spectral  properties  of  target  response 

The  analysis  till  now  has  concentrated  on  the  detection  of  a  target,  and  the 
determination  of  depth  and  principal  polarizabilities  of  the  target.  In  the  introduction  we 
also  described  the  vital  step  of  measuring  the  polarizabilities  as  a  function  of  frequency,  or 
time,  to  determine  the  shape,  size,  metal  content  and  even  the  wall  thickness  of  the  object. 


50 


Only  the  variation  of  the  induced  moments  with  frequency  permit  the  determination  of  the 
various  shapes  shown  in  the  cartoon  of  Figure  3.1.1.  This  analysis  immediately  raises  the 
issue  of  the  necessary  bandwidth.  To  get  an  idea  of  the  bandwidth  over  which  the  principal 
polarizabilities  undergo  their  defining  variation  we  again  turned  to  the  simulator  to  find  the 
spectral  and  transient  responses  for  some  typical  UXO. 

The  response  of  a  target  is  defined  here  as  the  secondary  field  (B)  at  a  given 
receiver  for  an  incident  field  from  a  given  transmitter.  The  response  is  thus  a  function  of 
the  T-R  pair  as  well  as  the  properties  of  the  target.  Induction  coil  sensors  typically  measure 
the  time  derivative  of  the  secondary  field  B  so  the  dB/dt  response  is  often  used. 


Figure  3.3.1:  Normalized  secondary  fields  (real  and  imaginary  components)  as  a  function 
of  frequency  for  37  mm  aluminum  spherical  shell  of  various  thicknesses  at  0.75  m 
depth. 

To  illustrate  the  spectral  response  of  a  variety  of  targets,  we  have  chosen  to  use  a 
simple  horizontal  loop  transmitter  with  an  in-loop  vertical  receiver  deployed  directly  above 
the  target.  This  is  basically  the  model  for  the  EM61  commercial  system.  This 
configuration,  the  target  size  and  shape,  and  the  separation  of  the  T-R  system  from  the 
target  are  shown  to  the  right  of  the  response  plots  in  the  following  figures.  In  all  cases  the 
secondary  fields  in  nano-Tesla  (nT)  are  normalized  by  the  transmitter  moment.  We  have 
plotted  the  frequency  response  real  (or  in-phase)  and  imaginary  (out  of  phase  or 
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quadrature)  and  the  transient  response  for  B  and  dB/dt  to  illustrate  diagnostic  behaviors  in 
both  domains.  We  have  studied  the  response  of  non-magnetic  conductors  as  well  as  the 
more  common  permeable  conductors  and  for  the  simple  sphere  we  have  analyzed  the  effect 
of  changing  shell  thickness  on  the  response. 

The  first  analysis  shows  the  effect  of  shell  thickness  on  the  frequency  response  of 
an  aluminum  37  mm  spherical  shell  0.75  m  below  the  T-R  pair.  For  any  body  the  high 
frequency  response,  so  called  inductive  limit  response,  depends  only  on  the  size.  In  the 
frequency  domain,  variations  in  the  shell  thickness  produce  characteristic  changes  in  real 
and  quadrature  response  at  frequencies  in  the  decade  above  and  below  the  frequency  of 
peak  quadrature  response.  In  Figure  3.3.1  this  band  is  between  30  and  3000  Hz.  The  high 
frequency  could  be  estimated  at  a  frequency  of  10  kHz  so  an  ideal  bandwidth  for 
identifying  this  target  would  be  30  Hz  to  10  kHz  or  2.5  decades. 

In  the  time  domain,  Figure  3.3.2,  the  step  function  response  in  B  clearly  resolves 
thicknesses  from  the  early  time  asymptote  at  3x1 0'5  sec  (which  also  determines  the  size)  to 
about  3x1 0'3  sec.  The  thinnest  shell  has  a  pure  exponential  decay  while  the  solid  sphere 
has  a  ‘stretched’  response  becoming  exponential  only  beyond  10~3  sec.  The  sensor  need 
only  have  a  dynamic  range  of  1.5  decades  to  resolve  the  response. 

The  time  derivative  dB/dt,  measured  by  standard  induction  coil  sensors,  is  shown  in 
Figure  3.3.3.  Diagnostic  changes  in  the  response  occur  between  3x10'  and  10'  sec  —  a 
much  wider  bandwidth  than  that  for  B  and  a  dynamic  range  of  at  lest  2.5  decades  is 
required  to  define  the  responses.  The  time  rate  of  change  dB/dt  can  be  measured  in  units  of 
Tesla/sec.  A  more  convenient  unit  of  dB/dt  is  the  Volt/m2.  It  can  be  shown  directly  from 
the  Maxwell  equations  that  1  Volt/m  =  1  Tesla/sec  exactly. 
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Figure  3.3.2:  Normalized  magnetic  field  response  as  a  function  of  time  for  37  mm 
aluminum  spherical  shell  of  various  thicknesses  at  0.75  m  depth. 
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Figure  3.3.3:  Normalized  dB/dt  response  as  a  function  of  time  for  37  mm  aluminum 
spherical  shell  of  various  thicknesses  at  0.75  m  depth. 


The  frequency  domain  response  for  the  same  spherical  shells  but  with  a  magnetic 
permeability  of  200  is  shown  in  Figure  3.3.4.  The  permeability  introduces  the  opposing 
static  magnetization  which  drives  the  real  response  negative  at  low  frequencies  and  causes 


small  perturbations  in  the  quadrature  response  also  at  low  frequencies.  Resolution  of  shell 
thickness  for  a  magnetic  target  is  only  possible  at  frequencies  below  about  300  Hz.  This 
simple  model  provides  the  first  evidence  that  it  may  be  difficult  if  not  impossible  to  resolve 
shell  thickness  in  magnetic  targets  in  the  frequency  domain. 

The  prospect  is  improved  considerably  in  the  time  domain  for  B  as  shown  in  Figure 
3.3.5.  As  expected,  the  thickness  variations  are  manifested  at  late  time,  between  10'5  sec 
and  10 3  sec.  The  thinnest  shell  develops  an  exponential  decay  by  3x1  O'5  sec  while  the 
solid  hasn’t  become  exponential  by  10"3  sec.  The  intervening  thicknesses  could  be  well 
resolved  with  a  dynamic  range  of  2  decades.  The  resolution  of  shell  thickness  with  dB/dt  is 
markedly  less  than  with  B,  Figure  3.3.6.  Curve  separation  is  clear  over  only  one  decade  of 
time,  10"3  to  10"2  sec,  and  over  3  decades  of  amplitude  but  the  diagnostic  separation  only 
begins  after  over  two  decades  of  amplitude  decay  from  shut-off. 
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Figure  3.3.4:  Normalized  magnetic  field  response  as  a  function  of  frequency  for  37  mm 
magnetic  spherical  shell  of  various  thicknesses  at  0.75  m  depth. 
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Figure  3.3.5:  Normalized  magnetic  field  response  as  a  function  of  time  for  37  mm 
magnetic  spherical  shell  of  various  thicknesses  at  0.75  m  depth. 


time  (s) 


. ►' 

37mm 


- 1=0.5  mm 

- 1=1  mm 

— 1=3  mm 

- 1=10  mm 

- solid 


Figure  3.3.6:  Normalized  dB/dt  response  as  a  function  of  time  for  37  mm  magnetic 
spherical  shell  of  various  thicknesses  at  0.75  m  depth. 


These  sample  models  illustrate  some  general  practical  conclusion:  B  rather  than 
dB/dt  is  more  diagnostic  of  shell  thickness  for  magnetic  and  non-magnetic  objects,  and 
shell  thickness  of  magnetic  objects  is  difficult  to  resolve  in  the  frequency  domain. 
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The  general  analysis  for  optimizing  the  T-R  configuration  in  Section  3.2  used 
idealized  dipoles  (spheres)  as  the  target.  It  is  important  to  know  how  well  these  results 
apply  to  the  case  of  actual,  usually  elongate,  UXO.  In  Figure  3.3.7  we  have  shown  the 
simulated  results  for  an  EM61  style  in-loop  system  over  a  typical  37  mm  shell  0.75  m 
below  the  system.  The  shell  is  solid,  steel,  and  has  an  aspect  ratio  (length/diameter)  of  3:1. 
The  plot  shows  the  dB/dt  response  for  both  the  horizontal  and  vertical  orientation  of  the 
shell.  For  comparison  the  response  of  the  37  mm  sphere  is  included.  It  is  immediately 
evident  that  the  actual  shell  response,  for  both  orientations,  is  larger  than  the  sphere 
responses.  It  may  be  concluded  that  all  the  detectability  analyses  of  Section  3.2  are  worse 
case  scenarios  and  so  are  excellent  design  guides  for  a  working  system. 

In  Figure  3.3.7,  and  subsequent  plots  for  other  targets,  we  have  also  plotted  the 
response  from  the  conductive  ground  in  which  the  target  is  immersed.  The  responses  for 
two  ground  resistivities,  10  and  100  Ohm-m,  are  plotted.  The  ground  response  basically 
imposes  an  early  time  limit  on  the  time  window  available  for  target  discrimination.  Once 
the  target  response  falls  below  the  ground  response  it  will  be  poorly  resolved,  especially 
since  the  ground  response  itself  will  variable  due  to  the  inhomogeneous  nature  of  the  near 
surface.  (The  role  of  magnetic  ground  is  still  being  analyzed  and  will  be  described  in  a 
later  report.)  For  a  conservative  design  approach  we  assume  a  ground  of  10  Ohm-m. 

37  mm  ‘shell’  3:1  aspect  ratio  at  0.75m  depth 
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Figure  3.3.7:  Amplitude  of  dB/dt  response  for  37  mm  sphere,  horizontal  and  vertical  shells 
3:1  aspect  ratio  at  the  depth  of  0.75  m  as  a  function  of  time  together  with  responses  for 
10  Q-m  and  100  Q-m  half-space. 
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For  this  EM61  simulation  we  can  also  include  the  normalized  practical  noise  level 
for  this  system  at  610  psec.  Since  the  noise  level  is  not  known  at  other  times  into  the 
transient,  we  have  simply  assumed  that  the  noise  is  constant  at  this  level  for  all  times  and 
plotted  this  horizontal  noise  floor  accordingly. 

The  early  time  limit  imposed  by  the  ground  and  the  later  time  limit  imposed  by  the 
noise  sets  the  time  window  over  which  the  response  can  be  measured.  For  this  target  the 
window  is  roughly  from  2x10'  to  3x10'  sec  (for  the  horizontal  shell).  Referring  back  to 
Figure  3.3.6,  it  can  be  seen  that  this  window  is  adequate  to  define  the  size  (early  time 
asymptote)  of  the  sphere  and  to  resolve  its  thickness,  but  the  dynamic  range  is  almost  five 
decades  above  the  noise  floor.  This  may  be  difficult  to  achieve  in  practice.  Using  only  the 
later  time  data,  say  from  2x1  O'4  sec,  may  be  satisfactory  especially  since  it  should  be  noted 
that  the  orientation  has  already  been  determined  in  the  first  stage  of  the  inversion. 

The  response  of  this  elongate  target  also  reveals  the  fact  that  the  step  response  is 
different  for  the  different  polarizabilities  of  a  non-spherical  target.  The  responses  of  the 
horizontal  and  vertical  orientations  of  this  37  mm  shell  actually  cross  at  a  few 
microseconds,  the  horizontal  shell  giving  a  longer  response  at  very  early  time  than  the 
vertical  shell.  This  is  an  excellent  illustration  of  the  fact  that  the  ratio  of  the  polarizabilities 
is  not  the  same  as  the  geometric  aspect  ratio  of  the  body.  As  can  be  seen  in  Figure  3.3.7  the 
ratio  of  the  horizontal  and  vertical  responses  at  later  time  is  10:1.  (An  interesting 
observation  is  that  the  ratio  of  the  horizontal  responses  for  the  target  to  that  of  a  sphere  of 
the  same  diameter  is  3:1.  This  appears  to  hold  true  for  all  the  elongated  targets  we 
considered.  The  result  holds  for  the  B  response  as  well.) 

For  B  rather  than  dB/dt  the  window  where  the  target  response  exceeds  the  ground 
response  widens  to  3xl0'6  sec  at  the  low  end,  Figure  3.3.8.  Since  there  are  no  existing 
systems  that  measure  B  we  have  no  estimate  of  when  the  late  time  response  meets  the  noise 
floor.  However,  the  horizontal  and  vertical  responses  still  differ  by  a  factor  10  and  most 
importantly  the  observed  separation  (including  very  early  time)  is  confined  to  less  than  3 
decades  of  amplitude  variation.  Instrumentally  this  is  much  more  manageable  than  the 
much  higher  dynamic  range  required  for  a  dB/dt  system. 
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37  mm  ‘shell’  3:1  aspect  ratio  at  0.75m  depth 


1.75m 


Figure  3.3.8:  Amplitude  of  magnetic  field  response  for  37  mm  sphere,  horizontal  and 
vertical  shells  3:1  aspect  ratio  at  the  depth  of  0.75  m  as  a  function  of  time  together 
with  responses  for  10  Q-m  and  100  Q-m  half-space. 


Figures  3.3.9  through  3.3.14  show  the  B  and  dB/dt  responses  for  three  other  targets 
representative  of  the  extremes  of  UXO  to  be  characterized  -  a  22  mm  shell  with  an  aspect 
ratio  of  8:1  (Figures  3.3.9  and  3.3.10)  at  a  depth  of  0.75  m,  a  105  mm  shell  with  an  aspect 
ratio  of  4:1  at  a  depth  of  2.65  m  (Figures  3.3.1 1  and  3.3.12),  and  a  larger  155  m  shell  with 
an  aspect  ratio  of  4.4:1  at  a  depth  of  4.55  m  (Figure  3.3.13  and  3.3.14).  In  all  these 
simulations  the  ground  response  limits  the  early  time  response  to  microseconds  in  B  and 
10’s  of  microseconds  in  dB/dt.  The  dynamic  range  requirement  is  reduced  in  B  and 
responses  to  at  least  0.01  sec  (10  msec)  are  required  to  clearly  resolve  the  decay 
characteristics.  For  the  large  target  at  greater  depth  the  response  must  be  obtained  to  at 
least  100  msec  at  which  point,  at  least  for  this  sample  EM6 1-like  system,  the  response  in 
dB/dt  is  below  the  noise  level  (see  Figure  3.3.13). 
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22  mm  ‘shell’  8:1  aspect  ratio  at  0.75m  depth 


Figure  3.3.9:  Amplitude  of  dB/dt  response  for  22  mm  sphere,  horizontal  and  vertical  shells 
8:1  aspect  ratio  at  the  depth  of  0.75  m  as  a  function  of  time  together  with  responses  for 
10  Q-m  and  100  Q-m  half-space. 


22  mm  ‘shell’  8:1  aspect  ratio  at  0.75m  depth 
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Figure  3.3.10:  Amplitude  of  magnetic  field  response  for  22  mm  sphere,  horizontal  and 
vertical  shells  8:1  aspect  ratio  at  the  depth  of  0.75  m  as  a  function  of  time  together 
with  responses  for  10  Q-m  and  100  Q-m  half-space. 
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105  mm  ‘shell’  4:1  aspect  ratio  at  2.65m  depth 


Figure  3.3.11:  Amplitude  of  dB/dt  response  for  105  mm  sphere,  horizontal  and  vertical 
shells  4: 1  aspect  ratio  at  the  depth  of  2.65  m  as  a  function  of  time  together  with 
responses  for  10  Q-m  and  100  Q-m  half-space. 


105  mm  ‘shell’  4:1  aspect  ratio  at  2.65m  depth 


Figure  3.3.12:  Amplitude  of  magnetic  field  response  for  105  mm  sphere,  horizontal  and 
vertical  shells  4: 1  aspect  ratio  at  the  depth  of  2.65  m  as  a  function  of  time  together 
with  responses  for  10  Q-m  and  100  Q-m  half-space. 
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-dBz/dt  /  Mtx  (Tesla/s  /  (Amp-mA2)) 


155  mm  ‘shell’  4.4:1  aspect  ratio  at  4.55m  depth 


Figure  3.3.13:  Amplitude  of  magnetic  field  response  for  155  mm  sphere,  horizontal  and 
vertical  shells  4.4: 1  aspect  ratio  at  the  depth  of  4.55  m  as  a  function  of  time  together 
with  responses  for  10  Q-m  and  100  Q-m  half-space. 


155  mm  ‘shell’  4.4:1  aspect  ratio  at  4.55m  depth 
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Figure  3.3.14:  Amplitude  of  magnetic  field  response  for  155  mm  sphere,  horizontal  and 
vertical  shells  4.4: 1  aspect  ratio  at  the  depth  of  4.55  m  as  a  function  of  time  together 
with  responses  for  10  Q-m  and  100  Q-m  half-space. 
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These  examples  have  been  chosen  to  illustrate  the  properties  of  a  particular  system,  in 
this  case  an  EM61  simulation,  in  resolving  elongate  objects  at  a  variety  of  depths  and  in 
showing  the  time  window  constraints  imposed  by  the  ground  response  and  the  system  noise 
levels.  The  results  of  Section  3.2,  for  example  Figure  3.2.5  shows  the  dramatic 
improvement  in  detectability  that  occurs  when  a  multi-element  T-R  configuration  is  used 
on  a  grid  of  stations.  Figure  3.2.1  for  example  shows  that  the  155  mm  target  at  a  depth  of 
4.55  m  would  not  be  detectable  with  the  simple  in-loop  system.  This  of  course  is  also  seen 
in  Figure  3.3.13  where  only  the  vertical  155  mm  elongated  shell  rises  above  the  noise  in  the 
6xl0'4  sec  to  7xl0'3  sec  time  window.  Figure  3.2.5  indicates  that  the  three  transmitter  - 
three  receiver  configuration,  occupying  many  locations,  could  easily  determine  the  depth 
and  polarizabilities  to  better  than  10%  for  this  155  mm  target.  A  combination  of 
detectability  graphs  such  as  Figure  3.2.5  or  3.2.1 1  and  specific  response  plots  such  as  those 
shown  in  the  suite  of  Figures  from  3.3.4  to  3.3.14  establish  the  bandwidth  and  moment 
requirements  for  an  actual  field  system. 

In  summary,  it  appears  that  for  step  function  excitation,  transients  from  10 
microseconds  to  100  milliseconds  are  to  be  detected  for  a  practical  range  of  UXO  and  a 
practical  range  of  depths. 


3.4  The  transient  response  of  spheres  and  spheroids 

Results  presented  hitherto  are  results  of  an  algorithm  for  calculating  the  response  of 
a  conducting,  permeable  sphere  (including  a  spherical  shell  of  arbitrary  thickness)  in  either 
a  uniform  field  or  the  field  from  a  finite  source.  Much  of  our  early  analysis  on  the  depth  of 
detection,  array  configurations  etc.  was  based  on  this  spherical  target.  UXOs  are  generally 
not  spherical  and  so  we  devoted  considerable  effort  to  adapting,  and  improving  numerical 
codes  to  model  prolate  and  oblate  spheroidal  objects.  We  now  have  reliable  codes  for 
modeling  solid,  conducting  and  permeable  spheroids  with  aspect  ratios  (diameter  to  length) 
from  1:4  to  4:1.  For  certain  polarizations  of  incident  field  the  results  are  accurate  to  an 
aspect  ratio  of  10:1. 
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We  simulated  the  target  response  of  various  transmitter-receiver  configurations, 
with  step  function  excitation,  over  these  spheroids  for  various  attitudes  of  the  bodies.  For 
comparison  we  plotted  the  transient  responses  of  spheres  of  the  same  radius  as  the 
equatorial  radius  of  the  spheroid  and  discovered  that  the  transients  were  all  of  the  same 
shape  and  were  simply  shifted  versions  of  the  sphere  response.  A  typical  result  of  this 
analysis  is  shown  in  Figure  3.4.1.  Here  we  have  plotted  the  vertical  field  transient 
response,  in  dB/dt  for  two  targets  located  directly  below  a  horizontal  loop  transmitter  at  the 
depth  of  0.75  m.  The  plot  shows  the  transient  response  of  a  37  mm  (equatorial  diameter) 
prolate  spheroid  with  an  aspect  ratio  of  3:1.  The  length  of  the  spheroid  is  thus  111  mm. 
The  response  is  shown  for  the  spheroid  oriented  vertically  and  horizontally.  For 
comparison  the  responses  of  a  37  mm  and  1 1 1  mm  sphere  are  also  plotted.  Finally  the  two 
sphere  responses  are  re-plotted  after  each  is  multiplied  by  a  constant;  0.246  for  the  111  mm 
sphere  and  2.71  for  the  37  mm  sphere.  The  shifted  sphere  responses  lie  almost  exactly  on 
the  spheroidal  responses  over  the  whole  transient  window.  The  spheres  and  spheroids 
preserve  their  conductivities  and  permeabilities  through  this  scaling. 


3:1  aspect  ratio  at  0.75m  depth 


Figure  3.4.1:  Equivalent  sphere  response. 
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Smith  and  Morrison  (2005a)  have  presented  a  detailed  analysis  of  this  empirical 
relationship  and  showed  that  the  scaling  factors  are  simple  functions  based  on  the  volumes, 
dc  magnetizations  and  high  frequency  limit  responses  of  the  spheres  and  spheroids.  The 
scaling  is  valid  for  highly  conducting  bodies  with  relative  permeabilities  above  50  and  with 
aspect  ratios  up  to  4  or  5  to  one  -  in  other  words  most  UXOs. 

The  impact  of  this  discovery  on  interpreting  the  apparent  principal  polarizabilities 
obtained  from  the  data  is  enormous.  The  numerical  codes  for  obtaining  the  response  of  the 
ellipsoidal  body  are  complex  and  are  time  consuming  even  on  a  large  computer.  Fitting  an 
ellipsoid  to  the  data  in  an  inversion  scheme  using  these  numerical  forward  models  is  even 
more  time  consuming  and  certainly  not  practical  for  real  time  processing  during  a  survey. 
On  the  other  hand  the  algorithm  for  the  sphere  response  is  a  closed  form  analytic 
expression,  which  is  very  fast  to  execute.  Since  the  scaling  factors  are  simple  functions  the 
dimensions  of  the  actual  spheroid  are  consequently  derived  almost  instantly  from  the  two 
principal  apparent  polarizabilities.  Thus  the  important  properties,  size  and  true  aspect  ratio 
can  be  determined  very  quickly  from  the  data.  Finally,  inversion  for  the  conductivity  and 
permeability  can  likewise  be  done  very  quickly  for  the  equivalent  spheres  as  described  in 
Section  3.8. 

3.5  The  effect  of  a  repetitive  transmitter  waveform 

The  ideal  transient  is  modified  by  the  practical  need  to  use  some  sort  of  repetitive 
waveform  rather  than  the  unrealizable  infinite  step  function.  The  usual  choice  is  a  square 
wave  such  as  that  shown  in  Figure  3.5.1.  We  have  chosen  to  illustrate  this  analysis  using  B 
rather  than  dB/dt  because  the  results  are  more  easily  seen  in  the  transient  responses  in  B 
than  in  the  more  complex  and  rapidly  decaying  dB/dt  responses.  The  conclusions  are  the 
same  in  either  domain.  The  transients  in  B  in  the  off-time  of  the  waveform  are  plotted  for 
the  step  response  and  for  four  progressively  shorter  periods  of  the  repetitive  waveform.  As 
the  waveform  shortens  the  dc  component  is  removed  and  the  low  frequency  content  of  the 
step  function  is  cut-off.  Further,  there  is  the  effect  that  the  transient  from  one  turn  off  has 
not  recovered  before  it  is  hit  by  the  next  turn-on.  Clearly  the  character  of  the  transient  is 
dramatically  changed  and  the  evident  stages  of  the  decay  discussed  for  the  step  function  in 
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Section  3.4  above  are  no  longer  as  evident.  Consequently  the  deduction  of  the  physical 
properties  from  the  transient  is  compromised.  From  Figure  3.5.1  it  is  clear  that  the  longest 
period  possible  would  be  best.  However  for  optimum  signal-to-noise  we  shall  see  that  the 
waveform  must  be  averaged,  or  stacked,  and  so  the  longer  the  waveform  the  longer  the 
total  time  that  the  system  must  remain  at  a  measurement  point,  and  so  the  survey  speed  is 
reduced.  The  gain  in  information  from  lengthening  the  waveform  is  perhaps  only  a 
theoretical  gain  since  for  a  practical  transmitter  moment  the  transient  at  late  time  will  have 
fallen  below  the  noise  floor. 

One  way  to  address  this  problem  is  to  simply  plot  the  initial  amplitude  of  the  B 
transient  for  some  simple  transmitter-receiver  configuration  over  spheres  of  various  sizes 
and  at  various  depths.  The  amplitudes  for  four  representative  spheres  at  four  different 
depths  are  plotted  as  a  function  of  repetition  period  in  Figure  3.5.2.  This  range  of  target 
size  and  depth  encompasses  the  range  required  for  useful  UXO  detection  and 
characterization.  The  amplitudes  reach  an  asymptote  with  increasing  period,  the  period  at 
which  this  occurs,  Tc,  varies  from  about  10  2  seconds  for  the  large  sphere  to  10  4  seconds 
for  the  small  sphere.  The  transition  to  the  asymptotic  amplitude  is  slow  -  the  values  drop 
less  than  order  of  magnitude  when  the  period  shrinks  by  two  orders  of  magnitude.  For  a 

—3  _o 

practical  system  the  period  could  lie  between  10  and  10  seconds  (repetitive  rate  of  1000 
Hz  and  100  Hz). 
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Figure  3.5.1:  Transient  response  of  a  repetitive  square  waveform. 
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Figure  3.5.2:  The  amplitudes  for  four  representative  spheres  at  four  different  depths 
plotted  as  a  function  of  repetition  period. 

We  addressed  this  problem  in  a  more  formal  matter  by  using  the  forward  model 
simulator,  now  expanded  to  include  the  response  of  the  band  limited  receiver  and  the 
square  wave  repetitive  waveform,  in  an  inverse  mode  to  determine  the  accuracy  of 
recovery  of  the  physical  properties  for  various  periods  of  the  waveform  (Section  3.8).  This 
process  basically  allows  us  to  find  the  shortest  period  for  which  the  properties  can  be 
recovered  with  the  best  accuracy. 

For  this  analysis  we  used  the  repetitive  waveform  of  Figure  3.5.3.  The  analysis 
included  the  repetition  rate,  1/T,  and  the  relative  pulse  length,  a/b.  The  current  pulse 
amplitude  was  varied  such  that  the  product  of  the  current  squared  and  the  pulse  duration 
was  kept  constant  -  in  other  words  the  energy  per  pulse  was  the  same. 
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Figure  3.5.3:  Repetitive  waveform. 
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The  receiver  was  critically  damped  and  we  used  two  resonant  frequencies  (fo),  3.2 
kHz  and  33  kHz.  The  target  was  a  160  mm  diameter  sphere.  The  data  were  inverted  to 
determine  log(r),  the  log  of  the  ratio  of  the  conductivity  to  relative  permeability,  log(a/pr), 
and  the  log  of  the  product  of  conductivity  and  permeability,  log(ap).  As  discussed  in 
Section  3.6  the  product  op  is  not  well  resolved  from  practical  transient  data  because  the 
transient  cannot  be  obtained  at  times  large  enough  for  the  exponential  decay  to  have 
occurred.  However  for  large  op  even  if  roughly  estimated,  the  radius,  r,  and  the  ratio  of 
o/p  are  well  determined.  In  other  words  the  determination  of  o/p  is  very  weakly  dependent 
on  the  value  of  op. 

The  results  are  displayed  in  Figures  3.5.4  and  3.5.5  as  contours  of  the  uncertainty  in 
the  estimates  of  log(r),  log(o/p),  and  log(op)  as  a  function  of  the  repetition  rate  in  Hertz 
(on  the  horizontal  axis)  and  relative  pulse  length  (on  the  vertical  axis).  The  range  of 
repetition  rates  was  based  on  rough  analysis  such  as  that  in  Figure  3.5.2. 

Figure  3.5.4a,  b,  and  c  show  results  for  a  resonant  frequency  of  3.2  kHz.  For  both 
log(r)  and  log(o/p)  there  is  a  broad  ‘valley’  of  lowest  uncertainty  running  diagonally  up  to 
the  right  for  both  parameters.  These  results  suggest  that  an  optimum  repetition  rate  is  in  the 
10  -  100  Hz  range  for  relative  pulse  lengths  from  0.03  to  0.1  respectively. 

The  practical  choice  is  guided  by  the  voltage  and  current  waveforms  for  the 
constant  energy  pulse.  The  short  pulse  (a/b  =  0.03)  would  have  to  have  a  very  high  current 
with  accompanying  high  voltages  (V  =  L  *  dl/dt)  on  the  switch. 

The  exercise  is  repeated  for  a  resonant  frequency  of  32  kHz  and  the  results  are 
plotted  in  Figure  3.5.5a,  b,  and  c.  Here  the  repetition  rate  range  has  been  extended  to  slow 


67 


the  rapid  deterioration  in  the  estimates  of  log(r)  and  log(a/p)  above  a  few  hundred  Hertz. 
The  optimum  range  is  not  dramatically  different  -  a  repetition  rate  of  100  Hz  and  a  relative 
pulse  length  of  0. 1  would  still  be  adequate. 


(a) 


2  10  100  200 
Repetition  Rate  (Hz) 
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Repet  i  t  i  on  R  at e  ( H  z ) 


(C) 

Figure  3.5.4:  Uncertainties  in  (a)  log(r),  (b)  log  (a/p),  and  (c)  log  (ap)  for  160  mm  sphere, 
resonant  frequency  fo=  3.2  kHz,  and  critically  damped  receiver. 
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Figure  3.5.5:  Uncertainties  in  (a)  log(r),  (b)  log  (a/p),  and  (c)  log  (ap)  for  160  mm  sphere, 
resonant  frequency  fo=  33  kHz,  and  critically  damped  receiver. 


3.6  Inversion  for  size,  conductivity,  and  permeability 

In  preliminary  analysis,  data  from  a  set  of  receivers  and  transmitters  is  reduced  to 
estimates  of  equivalent  dipole  polarizability  as  a  function  of  time  in  (estimated)  principal 
coordinates  of  an  object,  with  accompanying  error  estimates.  Given  the  good 
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approximation  of  spheroid  responses  by  appropriately  scaled  sphere  responses,  axial  and 
transverse  spheroid  responses  can  be  interpreted  in  terms  of  equivalent  spheres.  Dipole 
polarizability  as  a  function  of  time  is  inverted  for  sphere  radius,  r,  conductivity,  a,  and 
relative  permeability,  pr,  using  the  empirical  distribution  evolutionary  algorithm  of  Smith 
et  al.  (1994),  minimizing  the  weighted  squared  misfit  between  the  measured  dipole 
polarizability  as  a  function  of  time  and  that  of  a  sphere  illuminated  with  a  uniform 
incident  magnetic  fields  with  the  same  waveform  as  used  for  the  measured  data,  and 
filtered  through  the  same  receiver  transfer  function  as  the  measured  data.  The  empirical 
distribution  evolutionary  algorithm  is  a  population  based  algorithm:  initially  a  wide 
variety  of  spheres  are  modeled  with  randomly  selected  parameters  r,  a,  and  pr,  and  in  a 
sequence  of  successive  iterations,  the  space  of  parameters  is  sampled  increasingly 
densely  in  the  vicinity  of  the  better  points  found.  It  is  a  real  parameter  genetic  algorithm 
with  arithmetic  recombination,  a  high  degree  of  tournament  based  elitism,  with 
perturbations  based  on  differences  of  retained  population  members  and  discarded 
population  members.  More  details  are  given  in  Smith  and  Morrison  (2004). 

Once  a  set  of  sphere  parameters  r,  a,  and  pr  minimizing  the  weighted  data  misfit 
is  found,  parameter  uncertainties  are  estimated  using  a  local  analysis  based  on  the  partial 
derivatives  of  the  data  with  respect  to  the  parameters  (the  Jacobian)  evaluated  at  the 
minimum.  Letting  F  be  the  Jacobian  for  the  weighted  data  with  respect  to  the  logarithmic 
parameters  log(r),  log(c)  and  21og(pr),  the  covariance  matrix  of  log  parameter 
uncertainties  is  (ftf)-\  with  squared  log  parameter  uncertainties  on  its  diagonal  and  log 
parameter  error  covariances  on  its  off-diagonals.  In  general,  we  find  that  estimated  errors 
in  log(a)  and  log(pr)  are  highly  correlated,  and  errors  in  log(r)  essentially  independent. 
More  specifically,  we  have  found  the  eigenvectors  of  the  covariance  matrix  correspond 
very  closely  to  log  parameter  combinations  log(r),  log(a)  -  log(pr)  =  log(a/pr),  and  log(a) 
+  log(pr)  =  log(apr),  with  log(r)  best  determined,  log(a/pr)  moderately  determined  and 
log(apr)  generally  very  poorly  determined.  We  find  that  log(opr)  only  becomes 
reasonably  well  determined  as  the  length  of  data  used  in  inversion  approaches  or  exceeds 
the  fundamental  time  constant  of  the  sphere  aprpor2/d2  (with  d  «  1.4271  for  highly 
permeable  spheres,  d  «  n  for  non-magnetic  spheres). 
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Figure  3.6.1:  Squared  weighted  data  misfit  as  a  function  of  trial  sphere  conductivity  and 
relative  permeability  for  bipolar  data  with  an  80  Hz  repetition  rate,  a  0.4  duty  cycle 
square  pulse  and  a  32  k  radian/s  resonant  frequency  critically  damped  receiver  using 
data  from  0.47  ms  to  3.75  ms,  for  synthetic  data  from  a  160  mm  diameter  steel  sphere. 

A  graphic  example  of  the  insensitivity  of  data  to  the  product  ap,  is  shown  in 
Figure  3.6.1,  which  shows  squared  weighted  data  misfit  as  a  function  of  trial  sphere 
conductivity  and  relative  permeability  for  bipolar  data  taken  with  an  80  Hz  repetition 
rate,  a  0.4  duty  cycle  square  pulse  (2.5  ms  each  half  cycle),  and  a  32  k  radian/s  resonant 
frequency  (<»o)  critically  damped  receiver  using  data  from  when  primary  field  transients 
are  assumed  to  have  sufficiently  decayed  at  0.47  ms  (=1 5/coo)  to  the  beginning  of  the  next 

n 

half  cycle  at  3.75  ms,  for  synthetic  data  from  a  160  mm  diameter  steel  sphere,  (a  =  10  / 
Ohm-m,  pr  =  180  ).  In  this  figure,  for  each  point  plotted  trial  sphere  a  and  pr  have  been 
specified  and  the  sphere  radius  found  minimizing  the  squared  weighted  misfit.  The  misfit 
surface  forms  a  very  long  narrow  valley  through  the  true  a  and  pr  values.  Local  analysis 
shows  that  the  product  apr  is  essentially  undetermined,  corresponding  to  moving  along 
the  bottom  of  the  valley,  however,  for  ap,  <  4  x  1 06  /  Ohm-m  the  valley  starts  to  slope  up 
steeply.  This  corresponds  to  trial  spheres  with  time  constants  less  than  2  ms  starting  to  fit 
the  data  poorly,  so  although  the  apr  product  is  indeterminate,  the  data  could  still  place  a 
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floor  on  permissible  opr  products.  Moving  from  the  minimum  at  right  angles  to  the 
valley  corresponds  to  varying  the  ratio  a/pr,  the  misfit  rises  quickly  as  the  a/p,  is  much 
better  determined. 

3.7  Choice  of  receiver  bandwidth 

The  analysis  of  the  role  of  pulse  length  and  repetition  period  assumed  a  perfect 
receiver,  that  is  a  receiver  that  did  not  introduce  additional  distortion  to  the  transient 
waveform.  Unfortunately  realizable  sensors  for  either  B  or  dB/dt  do  not  possess  the 
essentially  infinite  bandwidth  required  to  measure  the  transient  without  distortion.  Coils 
inherently  have  a  distributed  capacitance  so  that  the  voltage  output  is  only  strictly 
proportional  to  dB/dt  below  some  too,  called  the  resonant  frequency  of  the  coil.  Such  a 
sensor  is  insensitive  to  dc  and  acts  as  a  filter  on  dB/dt  distorting  the  transient.  Some  high 
frequency  limit  is  desirable  in  any  sensor  to  reduce  the  broadband  ambient  electromagnetic 
field  noise.  For  practical  reasons  there  is  no  point  in  admitting  noise  of  frequencies  greater 
than  the  highest  frequency  needed  to  accurately  represent  the  transient  while  maintaining  a 
good  signal-to-noise  ratio.  Some  balance  must  be  found  between  moving  the  resonant 
frequency  to  a  high  value  to  accurately  represent  the  transient  while  letting  in  significant 
noise  to  moving  it  to  a  low  value  and  reducing  noise  but  significantly  distorting  the 
transient. 

A  variant  of  this  induction  coil  sensor  uses  a  feedback  winding  around  the  primary 
sensor  coil.  The  output  voltage  is  fed  back  to  this  auxiliary  winding  through  a  resistor  in 
such  a  way  as  to  create  a  field  through  the  primary  sensor  which  is  equal  and  opposite  to 
the  external  field  being  sensed.  The  current  in  the  feedback  winding  is  thus  largely 
proportional  to  the  sensed  field  rather  than  its  derivative.  These  feedback  sensors  have 
been  used  for  many  years  in  geophysical  exploration  systems  and  are  usually  described  as 
B  sensors.  The  feedback  is  not  ideal  and  the  resulting  sensor  response  is  bandwidth 
limited.  The  bandwidth  cannot  obviously  extend  to  dc  and  it  would  be  undesirable  to  let  in 
noise  above  those  frequencies  present  in  the  desired  transient.  Thus  the  B  sensor  also  has  a 
bandwidth  which  distorts  the  transient. 
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To  illustrate  this  discussion  we  consider  two  possible  receivers.  The  corresponding 
bandwidths  are  shown  schematically  in  Figure  3.7.1.  The  ideal  infinite  bandwidth  is 
represented  by  the  horizontal  line.  The  B-field  coil  (Figure  3.7.1a)  has  a  center  frequency, 
fo,  of  3  kHz  (this  is  actually  the  resonant  frequency  of  the  primary  sensing  coil)  and  has  flat 
B  response  from  fo/100  to  100*fo.  The  critically  damped  dB/dt  receiver  (Figure  3.7.1b)  has 
a  resonant  frequency  of  300  kHz  and  its  bandwidth  is  from  (v2  -1)  f0  to  (a/2  +1)  fr. 
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Figure  3.7.1:  Schematic  of  the  bandwidths  for  (a)  B  sensor,  and  (b)  critically  damped 
dB/dt  sensor. 


We  first  note  that  there  are  two  phenomena  that  have  a  major  impact  on  the 
observed  transient.  The  first  is  the  bandwidth  distortion  or  filtering  effect  described  above. 
There  is  also  the  fact  that  the  receiver  is  subject  to  the  sharp  transient  turn-on  and  turn-off 
of  the  transmitter.  The  receiver,  i.e.  the  physical  coil  of  wire,  develops  an  enormous  emf 
during  the  transmitter  transitions  which  persists  in  the  coil-amplifier  circuit  after  the 
cessation  of  the  transmitter  current.  In  a  critically  damped  dB/dt  receiver  this  transient 
signal  decays  with  the  time  constant  of  the  coil  -  equal  to  the  reciprocal  of  the  resonant 
frequency  of  the  coil,  i.e.  1/coo-  The  amplitude  of  this  “equipment”  transient  is  enormous  so 
it  is  likely  to  dominate  any  response  from  a  buried  target  for  a  considerable  time  after  the 
turn-on  or  turn-off  at  the  transmitter.  This  problem  could  be  avoided  by  using  a  receiver 
coil  oriented  orthogonal  to  the  primary  field  of  the  transmitter  -  the  so-called  minimum 
coupled  configuration.  In  practice  however,  this  is  difficult  to  do  because  the  orthogonality 
must  be  achieved  to  milliradians  or  less  and  maintained  at  this  level  as  the  platform  is 
moved  over  the  surface. 

Unfortunately  the  problem  is  exacerbated  in  the  B,  or  feedback,  receiver  where  the 
primary  field  transient  is  smaller  but  still  dominates  the  target  signal  by  many  orders  of 
magnitude  throughout  the  observation  interval. 
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The  following  brief  analysis  of  the  response  of  a  typical  induction  coil  sensor  to 
secondary  field  excitation  shows  in  greater  detail  the  role  of  the  detector  characteristics  in 
the  distortion  of  the  target  response. 

A  typical  TEM  system  is  shown  in  Figure  3.7.2.  It  consists  of  a  transmitter  with  a 
“step-off’  current  waveform,  an  induction  sensor  defined  by  its  electrical  parameters  L,  r, 
C,  and  a  “loop”  target  defined  by  its  time  constant  t.  We  will  assume  unity  gain  for  the 
sensor  amplifier  so  that  the  relationship  between  the  sensor  output  and  the  secondary 
magnetic  flux  produced  by  the  target  is  given  by: 

Vo  (t)  =  -dfdt  (3.7.1) 

For  the  purpose  of  the  following  discussion  we  will  assume  that  the  system  is  null  coupled 
so  that  4>o  (t)  =  0.  In  the  same  instances  it  is  convenient  to  normalize  the  sensor  output 
voltage  by  the  effective  area  of  the  sensing  coil,  Aef.  In  that  case,  since  <|)  =  BAef  we  get 
Vo/Aef  =  -dB/dt.  Once  again  we  see  that  1  nV/m2  =  1  nT/sec. 


Transmitter  Sensor 


Conductor 


Figure  3.7.2:  A  typical  TEM  system:  It  consists  of  a  transmitter  with  a  “step-off’  current 
waveform,  an  induction  sensor  defined  by  its  electrical  parameters  L,  r,  C,  and  a 
“loop”  target  defined  by  its  time  constant  i. 
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3.7.1  The  critically  damped  detector 


All  conventional  TEM  systems  use  critically  damped  induction  coil  sensor  to  detect 
the  target  signal.  It  is  either  self-resonant  or  tuned  to  a  central  frequency  ooo-  For  this  type 
of  sensor  shown  in  Figure  3.7.2  the  resonant  frequency,  fo,  is  given  by: 

2  TT  fo  =  coo  =  (LC)~1/2  (3.7.2) 

while  critical  damping  is  assured  by  having 

cd0  L/r  =  1/2  (3.7.3) 

In  that  case  the  detector  response  to  a  sinusoidal  flux  signal  §  =  <|)o  sin(cot)  is  given  by: 

V 

4>o®o 

As  shown  in  Figure  3.7.3  the  sensor  bandwidth  defined  by  the  half-power  points  equals  2(o0 
and  spans  the  range: 

(V2  -  1)  OOq  <  co  <  (V2  + 1)  oo0 

In  this  illustration  co0  is  set  to  105  s'1  corresponding  to  a  central  frequency  of  about  16  kHz 
so  that  the  sensor  has  a  bandwidth  of  about  32  kHz. 


CO/  COq 


(1  +  co/  conV 


(3.7.4) 


(Do  [S1] 

Figure  3.7.3:  Bandpass  characteristics  for  a  critically  damped  detector  from  Figure  3.7.2 
with  a  central  frequency  too  =  105  s'1  (»16  kHz). 
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The  time  domain  step  response,  i.e.  the  sensor  response  to  the  primary  flux  for  this 
detector  is  shown  in  Figure  3.7.4.  For  a  unit  step  in  the  ambient  flux  this  is  given  by: 

V  =  (<|>oa>o)  co0t  e  0)01 

It  rises  from  zero  to  a  maximum  value  of  10  psec  corresponding  to  a  detector  rise  time  of  5 
psec.  In  spite  of  this  fast  detector  response  it  should  be  noted  that  for  a  typical  maximum 
coupled  system  no  observations  of  the  secondary  field  can  be  made  until  the  primary  field 
“leakage”  i.e.  the  detector  response  falls  by  about  four  orders  of  magnitude  which  in  this 
case  occurs  at  about  100  psec. 


G>0  =  105  s1 


Figure  3.7.4:  Sensor  response  to  a  step  transition  in  primary  flux  for  the  induction  detector 
from  Figure  3.7.2. 

Now  let  us  turn  to  the  effects  of  bandwidth  on  signal  distortion.  This  is  the  only 
bandwidth  related  effect  seen  where  the  transmitter  and  receiver  are  in  minimum  coupling. 
Let  us  consider  a  target  signal  which  decays  with  a  time  constant  x: 

St  =  4>s  e_t/x  (3.7.5) 

In  this  case  the  ideal  induction  detector  would  produce  a  voltage  given  by: 
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Vj/(|)s  =--  e“t/T  +  8(t / x)  (3.7.6) 

x 

where  5(t/x)  is  a  delta  function  at  signal  onset  which  is  related  to  the  leading  edge  of  the 
transient  signal.  In  normalized  dimensionless  form  this  can  be  written  as: 

v  p_t^T 

V'  =  — ^-  = -  (3.7.7) 

<l>s®0  tX 

with  a  =  coot.  On  the  other  hand  the  voltage  output  for  a  practical  critically  damped 
detector  is  given  by: 


u  _t/T  ix  tx 

- Te  +  - - o 

(a-1)2  a-1  (a-1)2 


(3.7.8) 


The  first  term  of  this  expression  corresponds  to  the  decaying  signal  somewhat  amplified  by 


a  factor  of  g  = 


.  The  second  term  is  the  detector  response  to  the  leading  edge  of 


the  signal  (i.e.  its  initial  value).  It  is  a  distorted  version  of  the  delta  function  part  of  the 
ideal  detector  response  and  closely  resembles  the  actual  detector  step  response.  [The  step 
response  of  an  ideal  detector  is  a  delta  function.] 

The  “ideal”  and  “actual”  observed  signals  are  shown  in  Figure  3.7.5  for  a  signal 
time  constant  of  100  psec  and  a  critically  damped  detector  resonant  frequency  of  about  16 
kHz  (©o  =  105  s'1).  The  left  hand  side  shows  the  signal  distortion  phenomenon  in  true 
shape  on  linear  scales.  As  predicted  the  initial  part  of  the  observed  distorted  signal  indeed 
resembles  a  smeared  out  delta  function  while  the  late  part  of  the  observation  is  a  scaled 
(augmented)  version  of  the  true  signal.  As  shown  in  the  right  hand  portion  of  the 
illustration  the  signal  distortion  disappears  at  a  time  of  about  100  psec.  After  this  time  one 
can  extract  the  true  time  constant  of  the  signal  and,  if  desired,  correct  the  observations  by 
the  appropriate  “g”  factor  to  extract  the  correct  value  of  the  target  moment. 
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Figure  3.7.5:  Distortion  of  the  secondary  field  transient  for  a  100  (as  target  by  a  critically 
damped  detector  tuned  to  about  16  kHz  on  linear  and  on  log-linear  scales. 

We  will  call  the  instant  at  which  the  observed  and  augmented  target  signals 
coincide  to  within  a  given  percentage  the  early  observation  time  te.  It  can  be  calculated 
numerically  from  equation  (3.7.8).  The  results  for  our  sample  100  psec  target  are  shown  in 
Figure  3.7.6.  Here,  for  co0  =  105,  we  note  that  we  can  observe  the  augmented  transient  to  an 
accuracy  of  10%  at  about  80  psec  from  the  beginning  of  the  transient  while  if  an  accuracy 
of  1%  is  needed  then  the  earliest  observation  time  is  increased  to  1 10  psec.  One  can  derive 
a  general,  empirical  relationship  between  the  earliest  observation  time  for  the  undistorted 
transient  and  the  time  constant  (T=L/r)  of  the  critically  damped  detector.  For  10% 
accuracy  and  situations  such  as  this  one,  where  coq  ^  x ,  the  earliest  time  only  depends  on 

the  receiver  time  constant  (T  =  Vi  ooo)  and  is  approximately  given  by  te  =  20T.  Thus  the 
earliest  time  at  which  a  signal  can  be  correctly  observed  is  at  least  equal  to  twenty  times  the 
receiver  time  constant.  This  is  so  in  the  best  of  circumstances  where  the  receiver  time 
constant  is  already  twenty  times  smaller  than  that  of  the  target  and  we  have  a  minimum 
coupled  system. 


78 


Early  time 


Figure  3.7.6:  Earliest  time  of  detection,  to  1%  or  10%  accuracy  for  100  (as  target. 

We  thus  conclude  that  at  best  the  distortionless  observation  of  the  target  signal  is  a 
difficult  task. 

In  general  the  distortion  introduced  by  the  finite  bandwidth  of  the  sensor  cannot  be 
corrected  by  applying  a  system  correction.  The  data  outside  the  pass  band  is  reduced 
below  the  system  noise  and  deconvolution  operations  simply  bring  back  unacceptable 
noise.  Sensor  bandwidth  is  reduced  as  much  as  possible,  often  using  additional  filtering  in 
the  following  signal  processing  stages,  to  limit  the  noise  but  this  carries  the  problem  of 
distorting  the  transient  response.  Too  narrow  a  band  can  distort  the  signal  to  the  point 
where  it  is  actually  useless  for  detailed  characterization  of  the  target. 

The  feedback  coil  is  ideally  suited  for  optimizing  the  bandwidth.  Depending  on 
circumstances,  the  sensor  can  be  tuned  and  critically  damped  so  that  it  shows  the 
conventional  dB/dt  behavior.  On  the  other  hand,  with  the  judicious  use  of  feedback,  the 
sensor  can  be  made  to  have  a  flat  response  over  many  decades  of  frequency  so  that  its 
behavior  is  more  akin  to  that  of  magnetometer  that  measures  B.  In  either  case  the  dc  target 
response  cannot  be  recovered. 

The  idea  for  a  sensor  with  wideband  flat  frequency  response  originated  in  France 
nearly  forty  years  ago  (Glerc  and  Gilbert,  1964).  As  shown  below, 
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This  induction  device  consists  of  a  principal  winding  denoted  by  its  inductance  “L”  and 
resistance  “r”.  It  is  tuned  by  the  capacitance  “C”  to  a  central  frequency  “coo”  The  variable 
bandwidth  feature  is  incorporated  by  feeding  back  the  amplifier  (gain  A)  output  through  a 
current  limiting  resistor  R  to  an  auxiliary  winding,  ”1”,  which  is  inductively  coupled  to  the 
principal  winding  by  the  mutual  inductance  “m”.  The  receiver  response  is  a  function  of  the 
central  frequency  and  the  available  feedback  current.  It  is  the  latter  that  controls  the  system 
bandwidth  as  defined  by  the  ratio  of  the  highest  to  the  lowest  frequencies  in  the  pass  band, 
©h  /  ©l-  More  precisely, 

©H/©L  =  ©omA/2R  =  n2.  (3.7.9) 

The  frequency  response  for  a  variety  of  feedback  settings  as  indicated  by  the 
parameter  labels  which  correspond  to  values  of  coh/©l  is  shown  below  in  Figure  3.7.8. 


Figure  3.7.8:  Frequency  response  of  the  wideband  sensor. 
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Here  we  see  how  the  sensor  response  varies,  with  the  feedback  parameter  “n  ”,  as  a 
function  of  the  normalized  frequency  co/ooo-  A  feedback  parameter  value  of  unity 
corresponds  to  the  classical  critically  damped  induction  receiver  with  no  feedback.  If  cc>o, 
the  resonant  frequency  of  the  receiver,  is  much  greater  than  the  frequency  of  peak  target 
quadrature  response  then  the  sensor  will  in  fact  output  a  very  close  replica  of  the  time 
derivative  of  the  target  signal.  On  the  other  hand,  a  sensor  that  is  made  to  have  a  flat 
response  over  four  decades  of  frequency  by  the  application  of  massive  feedback  (n  = 
10,000),  will  in  fact  have  an  output  voltage  that  is  closely  proportional  to  the  actual 
secondary  magnetic  field  generated  by  the  target. 

The  discussion  to  this  point  has  considered  only  the  goal  of  recovery  the  B 
response.  It  is  considerably  easier  to  recover  the  true  dB/dt  response  by  moving  the 
resonant  frequency  of  the  coil  to  a  frequency  much  greater  than  the  frequency  of  peak 
target  quadrature  response  and  using  a  feedback  parameter  of  unity. 

For  a  maximum  coupled  receiver  the  response  is  often  dominated  by  the  transient 
induced  in  the  receiver  by  the  primary  filed  shut-off.  This  phenomenon  is  illustrated  in 
Figures  3.7.9  and  3.7.10.  Here  the  sensor  is  in  maximum  coupling  to  a  1-m  transmitter 
with  a  moment  of  180  Am2.  The  target  is  a  37  mm  iron  sphere  located  at  75  cm  below  the 
transmitter  plane.  For  the  critically  damped  sensor  the  dB/dt  response  of  the  target  alone, 
i.e.  with  an  ideal  infinite  bandwidth  receiver,  is  denoted  as  ‘unfiltered’.  The  target  as  seen 
by  the  critically  damped  receiver  of  Figure  3.7.1b  is  denoted  as  ‘filtered’.  The  primary 
field,  equipment,  transient,  denoted  as  ‘primary  step’,  is  over  five  orders  of  magnitude 
greater  than  the  signal  transient  at  times  earlier  than  10-5  seconds.  The  distorting  effect  on 
the  signal  transient  however,  becomes  insignificant  for  times  greater  than  25  psec. 

For  the  B  sensor  the  primary  step  response  is  given  by: 


V  _  ^0B  (-co y,t  _e-®Ht 

Aef  - 


(3.7.10) 


where  coH  =  n  cc>o,  ©l  =  oo/n,  with  n  defined  in  equation  (3.7.9),  V  is  the  measured  signal, 
and  Aef  is  the  effective  area  of  the  sensing  coil. 
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Even  though  the  signal  is  now  proportional  to  B  over  a  large  part  of  the  frequency  range  it 
still  has  to  be  measured  in  Volt/m  or  Tesla/sec.  As  shown  in  Figure  3.7.10  for  the  B 
sensor,  with  the  bandwidth  shown  in  Figure  3.7.1a,  the  primary  step  transient  is  about  four 
decades  of  magnitude  greater  than  the  signal  transient  and  it  has  essentially  no  decay  up  to 
10  seconds.  Any  practical  sensor  would  lack  the  dynamic  range  to  separate  these  two 
phenomena. 


Figure  3.7.9:  Critically  damped  dB/dt  sensor  response. 


Figure  3.7.10:  Wideband  B  sensor  response. 
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As  might  be  expected  the  signal-to-noise  ratio  for  any  receiver  is  a  function  of  its 
resonant  frequency  and  bandwidth.  We  can  see  this  by  again  taking  the  situation  illustrated 
in  Figure  3.7.9  and  examining  the  effect  of  the  resonant  frequency  for  the  critically  damped 
receiver  on  the  signal-to-noise  ratio  (SNR)  and  the  required  dynamic  range.  For  the 
critically  damped  300  kHz  receiver  in  question  and  ignoring  system  noise  we  find  an  SNR 
of  about  25.  This  estimate  is  based  on  a  signal  observation  at  about  12  psec  after  the 
extinction  of  the  primary  field  where  the  target  signal  has  amplitude  of  about  25  pV/nr. 
For  this  bandwidth  the  natural  observable  noise  is  about  1  pV/m2. 

Lowering  the  receiver  resonant  frequency  to  30  kHz  would  indeed  reduce  the  noise 
by  a  factor  of  five  to  about  200  nV/m  but  then  the  target  signal  would  not  be  seen  until  130 
psec  after  extinction  where  its  amplitude  is  only  900  nV/m  .  In  this  case  lowering  the 
receiver  resonant  frequency  by  a  factor  of  ten  would  result  in  an  SNR  reduction  by  a  factor 
of  about  five. 

This  is  only  an  illustrative  example  done  for  external  noise  measured  at  a  relatively 
quiet  site.  In  a  presence  of  large  noise  it  is  quite  likely  that  the  30  kHz  receiver  frequency 
would  be  optimal.  In  either  case  the  required  dynamic  range  would  be  about  160  db.  The 
dynamic  range  requirement  can  be  much  relaxed  by  using  a  half  sine  transmitter  with  a 
much  higher  moment. 
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4.  LAB  PROTOTYPE 


4.1  Introduction 

The  analyses  developed  in  this  project  are  applicable  to  AEM  systems  on  any  scale. 
We  have  found  in  this  process  that  systems  scale  roughly  with  depth  of  burial  and  size  of 
the  target.  Thus  the  dimensions  of  the  transmitter  control  the  field  strength  at  depth  and 
spacing  of  receivers  controls  the  accuracy  of  depth  estimates.  On  the  other  hand  the 
detection  of  small  near  surface  objects  depends  on  having  enough  receivers  deployed 
spatially  so  as  not  to  miss  their  spatially  small  response  pattern.  It  appears  that  given  the 
possibility  of  an  arbitrary  number  of  receivers  that  the  overall  system  dimensions  should  be 
set  by  the  size  of  the  transmitter  loop  and  the  practicality  of  moving  this  loop  over  the 
ground. 

We  have  decided  in  this  preliminary  design  to  choose  the  popular  cart  mounted 
configuration.  With  a  horizontal  loop  dimension  on  the  order  of  1.0  square  meter  this 
configuration  has  been  shown  to  be  maneuverable  and  to  be  capable  of  exciting  responses 
from  typical  UXO  to  depths  of  a  meter  or  two.  In  a  sense  we  have  chosen  a  conservative 
design  based  on  systems  which  have  a  demonstrated  track  record  of  detecting  metallic 
objects  to  the  depths  of  interest.  The  important  advances  of  the  proposed  design  are  that 
the  data  can  be  used  to  characterize  the  target  and  in  so  doing  classify  it  as  possible  UXO  or 
as  scrap,  and  that  the  target  can  be  located  and  identified  from  a  single  position  of  the 
multi-sensor  system.  The  following  sections  summarize  the  considerations  included  in  the 
design  process. 

4.2  Prototype  system 

We  assembled  a  bench  prototype  system  illustrated  schematically  in  Figure  4.2.1. 
This  system  was  not  rigorously  optimized  along  the  lines  suggested  in  Section  3,  rather  it 
was  assembled  with  readily  available  components  to  get  a  working  system  whose 
performance  could  be  used  to  confirm  modeling  results  and  to  be  a  starting  point  for  an 
optimized  system. 
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Figure  4.2.1:  Schematics  of  the  bench  test  acquisition  setup. 


The  main  transmitter  loop  is  the  one  used  in  the  Geonics  EM-61  system.  It  is 
roughly  one  meter  square,  has  approximately  32  turns.  We  used  a  half  sine  of  19  amperes 
peak  for  a  moment  of  about  600.  We  synthesized  results  from  an  orthogonal  three-coil 
system  by  simply  moving  the  one  coil  to  three  positions  in  the  jig  shown  in  Figure  4.2.2. 

We  used  eight  “pancake”  receiver  coils  mounted  within  the  footprint  of  the 
horizontal  transmitter  loop.  The  receivers  were  originally  designed  for  an  EM  scale  model 
study  and  are  not  designed  for  this  system.  Nonetheless  they  are  adequate  for  this  proof  of 
concept.  In  Figure  4.2.2  these  receivers  are  the  round  wheel-like  shapes.  The  black  centers 
are  tapered  plastic  plugs,  which,  when  screwed  down  on  precisely  located  threaded  posts 
on  the  horizontal  platform  assure  that  the  coils  are  accurately  relocated  if  they  need  to  be 
removed.  The  choice  of  location  of  the  eight  receivers  was  made  following  the 
optimization  procedure  of  Section  3.2. 


Figure  4.2.2:  Photo  of  the  prototype  system. 
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The  current  source  is  a  repetitive  bipolar  half  sine  waveform  the  parameters  of 
which  are  shown  in  Figure  4.2.3.  Again  we  used  a  half  sine  pulser  developed  in  a  previous 
project.  Based  on  the  discussion  of  Section  3.5  the  repetition  rate  is  good  but  the  relative 
pulse  length  is  probably  too  large. 


Figure  4.2.3:  Half  sine  waveform  of  prototype  system  -  repetition  rate  of  90  Hz,  relative 
pulse  length  is  0.18,  and  maximum  current  A=  19  Amperes. 

Each  receiver  in  this  prototype  is  connected  to  a  channel  of  a  multi-channel 
oscilloscope  by  a  simple  twisted  pair  cable.  The  trigger  for  the  oscilloscope  is  yet  another 
receiver  coil  mounted  essentially  on  the  transmitter  winding.  The  scope  is  thus  triggered 
on  the  trailing  sharp  edge  of  the  half  sine  derivative.  Recording  is  delayed  by  20  psec  to 
avoid  switching  transients  and  the  equipment  transients  discussed  in  Section  3.4. 

The  transient  is  digitized  with  a  sampling  interval  of  2  psec.  The  digitized  transient 
is  averaged,  or  stacked,  over  100  records  and  the  averaged  transient  is  sent  to  a  computer 
for  storage  and  processing. 

The  system  has  a  response  from  all  metallic  objects  within  a  radius  of  two  or  three 
dimensions  of  the  object.  Unfortunately,  this  means  that  there  will  be  spurious  response 
from  the  laboratory  equipment,  nearby  filing  cabinets  or  rebar  in  the  floor  etc.  We  decided 
to  deal  with  this  problem  by  recording  a  set  of  averaged  transients  with  no  test  target 
present  and  then  subtracting  this  transient  from  that  obtained  with  the  target  present.  All 
the  transients  obtained  in  the  following  test  were  treated  in  this  manner. 

The  transmitter-receiver  system  was  mounted  about  one  meter  off  the  floor  and 
various  target  objects  were  placed  beneath  it.  The  individual  points  of  this  transient  were 
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further  averaged  into  time  windows  whose  width  is  equal  to  25%  of  the  center  of  the  center 
frequency  of  the  window.  This  reduced  the  number  of  points  representing  the  transient  to 
about  30.  The  averaged  transients  for  each  receiver,  and  for  three  orthogonal  positions  of 
the  transmitter  were  then  used  to  invert  for  the  location  and  principal  polarizabilities  (and 
the  orientation,  azimuth,  and  dip  of  the  major  polarizability). 

The  inversion  stopped  at  this  stage.  With  the  polarizabilities  determined  as  a 
function  of  time  it  is  possible  to  utilize  the  results  of  Sections  3.3  and  3.6  to  determine  the 
physical  aspect  ratio,  log(r)  and  log(a/p).  These  results  were  for  a  repetitive  square  wave 
and  algorithms  have  not  yet  been  modified  for  the  half  sine  waveform  used  in  the 
prototype. 

4.3  Test  model  results 

We  used  three  simple  targets  in  the  initial  tests;  a  simple  multi-turn  coil  18  cm  in 
diameter,  a  copper  sphere  138  mm  in  diameter  and  a  steel  ellipsoid  (Figure  4.3.1)  9  inches 
long  and  3  inches  in  diameter  (aspect  ratio  of  3: 1).  For  each  target  the  location  (x,  y,  z),  the 
three  principal  moments  and  the  orientation  of  the  three  principal  axes  are  plotted  as  a 
function  of  time  on  the  transient.  The  estimates  are  plotted  as  solid,  colored,  lines  and 
dashed  lines  of  the  same  color  indicate  the  variability  of  the  estimate.  On  each  plot  the 
actual  value  of  the  parameters  plotted  is  listed. 


Figure  4.3.1:  The  9  inch  (23  cm)  long  and  3  inches  (7.6  cm)  in  diameter  steel  ellipsoid. 

Figures  4.3.2a,  b,  c  are  plots  of  the  location,  principal  moments  and  principal  axes 
directions  (sometimes  given  as  direction  cosines)  for  the  18  cm  diameter  shorted  multi- turn 
loop.  The  loop  was  positioned  about  30  cm  below,  and  coaxial  with,  the  horizontal 
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transmitter  loop.  The  location,  Figure  4.3.2a,  is  recovered  very  accurately  between  5x10  5 
sec  and  5x10  4  sec.  The  principal  moment,  Figure  4.3.2b,  there  should  be  only  one  for  a 
planar  loop,  is  well  recovered  over  the  same  time  interval,  and  the  minor  moments,  which 
should  be  zero,  are  about  two  orders  of  magnitude  less  than  the  principal  moment  -  they 
are  in  the  noise  of  the  estimate.  Finally,  the  direction  cosines,  Figure  4.3.2c,  are  well 
recovered. 


Figure  4.3.2a:  Inversion  results  for  the  position  of  18  cm  diameter  shorted  multi-tum  loop. 


Figure  4.3.2b:  Inversion  results  for  the  principal  moments  of  18  cm  diameter  shorted 
multi-tum  loop. 
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Figure  4.3.2c:  Inversion  results  for  the  moment  orientations  of  18  cm  diameter  shorted 
multi-turn  loop. 

The  results  for  the  138  mm  diameter  copper  sphere  are  plotted  in  Figures  4.3.3a  and 
4.3.3b.  The  sphere  center  in  this  case  is  only  20  cm  (200  mm)  below  the  plane  of  the 
horizontal  coil  so  it  might  be  wondered  whether  the  dipole  representation  is  accurate. 
Indeed  the  inverted  parameters  at  ‘early’  time,  from  5x1 0~5  to  10~4  seconds,  seem  to  depart 
in  a  well  defined  way  from  the  true  values.  Nonetheless  the  location  is  well  recovered 
(Figure  4.3.3a),  and  the  principal  dipole  moments,  which  should  all  be  equal,  fall  roughly 
within  each  others  error  bounds  between  5xl0~4  and  10~4  seconds  (Figure  4.3.3b).  The 
moment  orientations  are  meaningless  in  this  case. 
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Figure  4.3.3a:  Inversion  results  for  the  position  of  138  mm  diameter  copper  sphere. 


Figure  4.3.3b:  Inversion  results  for  the  principal  moments  of  138  mm  diameter  copper 
sphere. 

Finally,  we  positioned  the  9  inch  (23  cm)  steel  ellipsoid  first  directly  below  the 
center  of  the  transmitter  at  the  depth  of  30  cm,  with  its  long  axis  horizontal  (dip  =  0°)  and  at 
an  angle  of  45°  to  the  principal  coordinates  of  the  system  (azimuth).  This  test  was  done  in 
an  all-wood  building  at  our  Richmond  Field  Station  and  in  general  the  noise  was  less.  The 
location,  Figure  4.3.4a,  is  recovered  very  well  over  the  entire  transient  window,  the 
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principal  polarizabilities  are  equally  well  resolved,  Figure  4.3.4b.  Finally  the  dip  and 
azimuth  are  recovered  very  well,  Figure  4.3.4c,  -  so  well  in  fact  that  it  appears  that 
ellipsoid  was  in  fact  2°  off  its  intended  45°  azimuth. 


Figure  4.3.4a:  Inversion  results  for  the  position  of  9  inch  (23  cm)  steel  ellipsoid. 


Figure  4.3.4b:  Inversion  results  for  the  principal  polarizabilities  of  9  inch  (23  cm)  steel 
ellipsoid. 
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Figure  4.3.4c:  Inversion  results  for  the  orientation  of  9  inch  (23  cm)  steel  ellipsoid. 

Although  the  system  is  by  no  means  optimized  and  the  noise  levels  are  well  above 
those  anticipated  in  a  final  system,  we  then  placed  the  ellipsoid  at  a  depth  of  45  cm  and 
again  horizontal  and  at  an  azimuth  of  45°.  This  time  the  center  of  the  ellipsoid  was 
displaced  25  cm  in  x  and  y  from  the  centerline,  or  z-axis.  These  positions  were  not 
measured  with  high  accuracy  so  the  location,  Figure  4.3.5a,  at  early  time  can  be  assumed  to 
be  more  accurate  than  the  intended  value.  The  recovered  polarizabilities,  Figure  4.3.5b, 
clearly  show  the  correct  apparent  aspect  ratio  but  there  is  a  difference  in  the  two  minor 
polarizabilities  that  is  larger  than  their  respective  errors.  We  do  not  understand  this  result 
at  present.  Early  time  estimates  of  the  dip  (error  less  than  ±1°)  and  azimuth  (error  less  than 
±5°),  Figure  4.3.5c,  are  excellent.  All  the  results  become  very  noisy  beyond  about  2x  10~4 
seconds  as  the  transient  signals  from  such  a  relatively  deep  target  drop  into  the  noise. 

These  results,  especially  the  last  for  an  off-center  oriented  ellipsoid  at  the  depth  of 
45  cm,  clearly  show  that  the  multi  element  system  can  detect  the  object  and  more 
importantly,  determine  its  principal  polarizabilities  and  their  directions.  Moreover,  it  is 
evident  that  these  quantities  can  be  obtained  over  a  time  window  that  will  allow  subsequent 
inversion  for  the  true  physical  aspect  ratio,  the  size  and  the  ratio  of  a/p. 


92 


Figure  4.3.5a:  Inversion  results  for  the  position  of  9  inch  (23  cm)  steel  ellipsoid  with 
offset. 


Figure  4.3.5b:  Inversion  results  for  the  principal  polarizabilities  of  9  inch  (23  cm)  steel 
ellipsoid  with  offset. 


Figure  4.3.5c:  Inversion  results  for  the  orientation  of  9  inch  (23  cm)  steel  ellipsoid  with 
offset. 
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5.  CONCLUSIONS  AND  ACCOMPLISHMENTS 


The  most  significant  finding  are  detailed  in  Section  2  of  this  report  where  we  use 
the  principal  dipole  moments  model  to  evaluate  the  improvements  in  system  performance 
that  can  be  obtained  with  use  of  multiple  receivers  and/or  transmitters.  The  other  most 
significant  result  here  is  that  the  depth  of  detection  can  be  doubled  compared  to  a 
conventional  coaxial  transmitter  -  receiver  system  if  one  uses  a  3  component  transmitter  as 
well  as  a  3  component  receiver. 

An  optimal  system  is  one  that  allows  the  most  accurate  location  and  classification 
of  the  target  at  the  lowest  cost  in  terms  of  the  system  size,  weight,  and  power.  Basically  we 
are  taking  the  “most  bang  for  your  buck”  approach.  Although  we  are  beginning  to  see 
traditional  optimization  possibilities  in  coil  size  and  system  bandwidth,  i.e.  a  particular 
setting  is  better  than  all  others,  a  more  limited  approach  may  have  to  be  taken  to  the  overall 
system  design.  This  implies  increasing  system  complexity,  weight,  and  power  to  the  point 
of  diminishing  returns  because  an  absolutely  optimum  value  for  some  of  the  parameters 
may  not  exist.  Thus  increasing  the  number  of  sensors  and/or  transmitters  beyond  a  certain 
level  may  not  detract  from  system  performance  but  does  not  improve  it  either. 

A  powerful  simulator  has  been  implemented  for  determining  the  optimum 
transmitter-receiver  configuration  for  UXO  detection  and  characterization.  Another 
simulator  has  been  developed  for  analyzing  the  role  of  bandwidth  and  noise  in  determining 
the  spectral  response  of  the  principal  dipole  moments. 

The  major  conclusions  from  research  and  development  are: 

1)  A  multicomponent  transmitter-receiver  system  is  essential  for  the  identification  of  the 
principal  dipole  moments  of  a  target,  be  it  UXO  or  clutter.  A  complex  target  can  be 
detected,  located,  and  defined. 

2)  A  three  component  loop  transmitter  with  5  small  vector  sensors  suffices  to  uniquely 
determine  the  principal  dipole  moments  of  a  target.  Stand-alone  equipment  recovers  all 
target  parameters  from  a  single  position. 

3)  Similar  detection/characterization  results  can  be  obtained,  with  different  configurations, 
using  deployment  on  a  2D  grid,  on  a  profile  line,  or  in  a  stand  alone  mode. 

4)  The  utility  of  grid  or  line  deployment  is  limited  by  positioning  accuracy. 
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5)  For  maximum  utility  a  system  capable  of  stand  alone  detection/characterization  seems 
best  because: 

i)  it  can  be  used  effectively  when  access  does  not  permit  precise  line  or  grid  operation. 

ii)  its  sensitivity  is  independent  of  absolute  position. 

iii)  when  used  in  line  or  grid  mode  its  sensitivity  is  superior  to  any  other  configuration. 

6)  )  A  dB/dt  sensor  is  preferable. 

7)  The  bandwidth  of  the  receiver  and  its  associated  electronics  have  a  profound  effect  on 
the  secondary  field  transient  and,  in-tum,  on  the  accurate  recovery  of  the  distinctive 
properties  of  the  transient  that  enable  the  identification  of  the  spectral  properties  of  the 
target. 

8)  A  bandwidth  of  at  least  four  decades  appears  to  be  necessary  to  describe  adequately  the 
spectral  response  of  the  principal  moments.  Low  frequencies  are  required  to  identify  the 
permeability  of  the  target  and  this  in-tum  requires  a  low  pulse  repetition  rate  for  the 
transmitted  fields. 

9)  The  large  bandwidth  admits  more  noise  and  thus  requires  higher  transmitter  power  to 
maintain  an  adequate  signal  to  noise  ratio. 

10)  Target  size,  ratio  o/pr,  and  tme  physical  aspect  ratio  can  be  determined. 

11)  The  ground  response  imposes  an  early  time  limit  on  the  time  window  for  target 
discrimination. 
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Figure  A.l.1.1:  Largest  26  eigenvalues  of  K  scaled  to  inverse  decay  times  X  as  a  function 
of  aspect  ratio,  for  prolate  spheroids  of  conductivity  a  =  3.5  10  Q  m'  and  with  a  6 
cm  equatorial  radius,  with  equation  (A.1.10)  truncated  at  polynomial  order  p=  7....  108 
Figure  A.l.1.2:  Receiver  coil  voltage  as  a  function  of  time  for  a  simulated  Geonics  EM-61 
system  with  coils  0.6  m  above  a  6  cm  radius  aluminum  sphere.  Previous  results  (large 
circles).  Modified  code  results  summed  over  196  modes  (dashed).  Analytic  (solid). 

. 109 

Figure  A.l.2.1:  EM-61  simulation . 110 

Figure  A.l.4.1:  Transient  responses  for  B  and  dB/dt  for  6  cm  conductive  permeable  steel 
sphere,  at  1  m  depth . Ill 


104 


Figure  A.l.4.2:  Transient  responses  for  B  and  dB/dt  for  6  cm  conductive  permeable  steel 


sphere,  for  varying  shell  thickness,  at  1  m  depth . 112 

Figure  A.  1.4.3:  Graphical  user  interface  for  calculating  responses  of  a  spherical  target 

with  an  arbitrary  transmitter-receiver  configuration . 1 13 
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APPENDIX  1:  NUMERICAL  MODELS 


A.1.1  Improvements  to  the  Weichman-Blackhawk  MFT  Code 

P.B.  Weichman  of  Blackhawk  Geophysics  made  available  source  code  for  an  integral 
equation  method  for  computing  the  secondary  (scattered)  fields  due  to  electromagnetic 
induction  in  conductive  and  possibly  magnetic  spheroids.  The  code  expands  Cartesian 
components  of  the  electric  field  (and  of  the  magnetic  field  when  necessary)  in  terms  of  a  set 
of  polynomials  f(r); 

oo  3 

E(r)  =  ZZrlij  fi(r)Xj  ’  (A.1.1) 

i=l  j=l 


where  Xj  is  a  unit  vector  in  the  j'th  Cartesian  coordinate  direction.  In  its  course  of 
computations,  the  code  finds  eigenvalues  and  eigenfunctions  of  an  equation  of  form 


E(r)  =  (a'-aj,) 
(in  Gaussian  units),  where 


4n  pbA.  VV  • 


j  gb(r,r')E(r)  dr' 

spheroid 


(A.  1.2) 


a  =  a- 


471 


(A.  1.3) 


and  gb(r,r')  is  the  Green’s  function  for  a  scalar  Helmholz  equation  in  the  background 
medium  (  Ob,  Sb,  pb).  Terms  involving  (p  -  pb)  H  (r' ),  arising  for  magnetic  spheroids,  have 
been  omitted  for  simplicity  of  exposition.  Integrals  of  form 

Ii(rk)  =  {  gb(rk>r')fi(r')dr'  (A.  1.4) 

spheroid 

are  evaluated  over  a  set  of  m  points  rk,  (k  =  l,m),  with  the  Green’s  function  approximated 
within  the  spheroid  as  l/(47i|r  -  r'|) ; 

w-wn-A  j  fif’) 


dr’  . 


(A.  1.5) 

y  —  ri 

""spheroid  I  I 

For  f  (r')  a  polynomial  of  order  p  in  x,  y,  and  z  (e.g.,  xa  y|!  zp'a'(3),  I';  (r)  is  a  polynomial  of 
order  p+2  or  lower  within  the  spheroid,  so  can  be  expanded  as 

V2 


i'=l 


(A.  1.6) 
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where  np+2  is  the  number  of  linearly  independent  polynomials  in  x,  y,  and  z,  of  order  p+2  or 
less.  Coefficients  yi'(l)  are  found  by  numerically  integrating  (A.  1.5)  at  the  points  rk,  (k  = 
l,m),  for  m  >  np+2,  and  inverting  (A.  1.6)  for  y,'*0,  i'  =  1,  np+2.  With  y/°  computed,  it  is 
straightforward  to  evaluate 

np+2 

VV-ij  Ij(r)  =  X  y!'0  VV-Xj  fj'(r)  .  (A.  1.7) 

i'=l 

This  is  re-expanded  in  the  form 


W-Xj  I-(r)  =  X  X7?jkfi'(r)xk  , 

i'=l  k=l 


(A.  1.8) 


by  expressing  the  needed  derivatives  of  each  function  f;(r)  in  terms  of  the  basis  functions 
fi(r).  Using  expansions  (A.  1 . 1),  (A.  1 .6),  and  (A.  1 .8)  in  equation  (A.  1 .2)  yields  (A.  1 .9) 

oo  3 


oo  3  oo  3 

ZZTiijfi(r)ij = (a'-  ab)  zz 

i=l j=3  i=l j=l L 


fi'Wij  -Vi  Er?jL  fr(r)i 

abi'=lk=l 


i'=l 


Truncating  the  expansion  to  the  np  polynomials  of  order  p  or  less  yields  (A.  1.10) 

nP  3 


nP  3  nP  3 

ZZTHifi(r)ij  =  (<7'-<7b)ZZ 

i=l j=3  i=l j=l 


^ErS^'Wij-VEEr^fiV) 

c  i'=l  ab  i'=i  k=l 


Pij  • 
(A.  1.9) 


hij  • 


(A.  1.10) 

As  originally  programmed,  integral  (A.  1.6)  was  fitted  by  least  squares  to  the  np 
polynomials  of  order  p  or  less,  at  points  rk,  k  =  l,m,  yielding  approximate  values  for  the 
first  np  coefficients  y,'(l).  As  a  consequence,  the  final  summation  over  i'  in  equation 
(A.  1.10),  was  missing  terms  of  orders  p-1  and  p  arising  from  the  VV-  term  in  equation 
(A.  1.2). 

Equation  (A.  1.10)  is  satisfied  throughout  the  spheroid  when  satisfied  separately  for 
the  coefficients  of  each  pair  f(r)  Xj ,  i  =  1,  np,  j  =  1,3,  yielding  a  matrix  equation  in  a  vector 


e  of  coefficients  q,, ; 

e  =  Ke  ,  (A.  1.11) 

(sic).  When  the  set  of  polynomials,  f(r),  is  orthogonal  over  the  spheroid  (and  normalized), 
the  resulting  matrix  K  is  symmetric,  so  has  real  eigenvalues  and  eigenvectors. 
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As  originally  programmed,  the  expansion  was  in  terms  of  monomials 
fj(r)  =  xa‘y^zYi  with  oq  +  P;  +  y;  <  p.  In  that  basis,  truncation  of  equation  (A.  1.9)  in 

equation  (A.  1.10)  and  formation  of  equation  (A.  1.11),  results  in  a  matrix  K  which  is  not 
symmetric,  and  has  complex  eigenvalues.  Keeping  terms  up  to  i  =  np+2  in  expansion 
(A.  1.6),  and  the  corresponding  terms  in  the  first  summation  over  i'  on  the  right  hand  side  of 
equation  (A.  1.10)  by  summing  to  i'=  np+2,  expressing  (A.  1.10)  in  a  basis  of  polynomials 
orthogonal  over  the  spheroid,  truncating  in  that  basis  to  polynomials  of  order  less  than  or 
equal  to  np,  and  re-expressing  in  the  monomial  basis,  is  sufficient  to  ensure  that  the 

eigenvalues  of  K  are  real  within  round-off  error.  It  is,  however,  simpler  to  leave  the 
problem  expressed  in  terms  of  the  orthogonal  basis,  as  then  subroutines  taking  advantage  of 

the  symmetry  of  K  can  be  used  to  find  its  eigenvalues  and  eigenvectors.  With  these 
changes,  the  leading  eigenvalues  are  stable  up  to  aspect  ratios  of  28:1  for  non-magnetic 
prolate  spheroids  (Figure  A.  1 . 1 . 1 ). 


Lowest  26  modes  for  Ellipsoid  with  6  cm  minor  axes 
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Figure  A.l.1.1:  Largest  26  eigenvalues  of  K  scaled  to  inverse  decay  times  A  as  a  function 
of  aspect  ratio,  for  prolate  spheroids  of  conductivity  a  =  3.5  107  Lf'm'1  and  with  a  6 
cm  equatorial  radius,  with  equation  (A.  1.10)  truncated  at  polynomial  order  p=  7. 
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Once  values  of  inverse  decay  time  X  are  found  for  which  matrix  K  has  unit 
eigenvalue(s),  the  corresponding  eigenvectors  give  the  expansion  coefficients  for  the 
corresponding  mode  of  the  spheroid.  The  excitation  of  the  modes  by  a  transmitter  loop  of  a 
specified  geometry  is  computed  for  a  ramp-on/ramp-off  transmitter  current,  and  the 
contributions  of  the  different  modes  to  voltages  observed  in  a  receiver  coil  are  summed 
over  modes.  Figure  A.  1.1. 2  shows  results  from  the  modified  code  (dashed  line)  using  the 
lowest  196  modes  for  a  6  cm  radius  sphere  of  conductivity  a  =  3.5  107  Q  'm"1  (with 
fundamental  time  constant  x«  15  ms)  at  12  cm  depth  below  a  simulated  Geonics  EM-61 
system.  Also  shown,  are  analytic  results  for  the  same  sphere  (solid  line),  and  results  from 
the  unmodified  code  (large  circles),  both  for  the  same  simulated  transmitter-receiver 
system.  The  modified  code  results  agree  to  within  10  %  0.5  ms  (x/30)  after  transmitter 
pulse  end,  and  to  within  2  %  3  ms  (x/5)  after  pulse  end.  In  contrast  results  for  the 
unmodified  code  are  accurate  to  within  20  %  4  ms  («  x/4)  after  pulse  end,  and  worse  before 
that. 


6  cm  radius  aluminum  sphere  60  cm  below  EM61 


Figure  A.l.1.2:  Receiver  coil  voltage  as  a  function  of  time  for  a  simulated  Geonics  EM-61 
system  with  coils  0.6  m  above  a  6  cm  radius  aluminum  sphere.  Previous  results  (large 
circles).  Modified  code  results  summed  over  196  modes  (dashed).  Analytic  (solid). 
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Using  the  modified  code  to  simulate  EM-61  responses  on  a  grid  of  system  locations 
above  a  8  cm  by  24  cm  aluminum  prolate  spheroid,  and  inverting  for  equivalent  dipole 
polarizability,  yields  an  equivalent  dipole  polarizability  that  is  consistent  with  the 
symmetry  properties  of  the  spheroid.  The  resultant  equivalent  dipole  polarizability  is 
centered  at  the  spheroid  center,  has  two  equal  principal  moments  and  a  third,  smaller, 
principal  moment  with  principal  direction  parallel  to  the  spheroid  symmetry  axis  (see 
Estimating  Equivalent  Dipole  Polarizabilities  Section  for  details). 

A.  1.2  Loop-Loop  over  Dipole  Targets 

We  have  developed  a  general  code  for  modeling  the  response  of  any  configuration 
of  finite  loops  over  a  target  represented  by  orthogonal  dipole  moments  of  any  magnitude. 
An  example  of  the  EM-61  response  at  early  time  over  a  simple  vertical  moment  target  (e.g. 
a  small  horizontal  circular  sheet)  is  shown  in  Figure  A.  1.2.1.  This  is  the  general  code  for 
modeling  specific  T-R  configurations  such  as  those  described  in  Section  3. 


x  (m) 

Figure  A.l.2.1:  EM-61  simulation. 

A.1.3  Sheet  in  layered  half  space 

We  have  included  a  finite  loop  source  in  our  general  code  for  determining  the 
response  of  a  finite  rectangular  thin  sheet  in  a  layered  and  conductive  half  space.  This  is  a 
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useful  general  model  for  evaluating  the  ratio  of  the  target  response  to  the  half  space 
response  and  we  have  used  it  to  show  that  this  ratio  decreases  as  the  transmitter  receiver 
separation,  or  the  size  of  the  transmitter  loop,  increases. 

A.1.4  Conductive,  permeable  spherical  shell 

We  have  completed  an  analysis  of  the  spectral  response  of  a  conductive,  permeable, 
shell  in  a  conductive  whole  space  using  the  closed  form  solution  we  programmed  in  the 
Phase  1  part  of  this  project.  The  source  is  either  a  uniform  field  (approximately  valid  when 
the  size  of  the  sphere  is  small  compared  to  the  distance  to  the  source)  or  the  field  of  a 
dipole  (and  by  superposition,  the  field  of  a  finite  loop). 

This  code  is  particularly  useful  in  showing  that  for  typical  UXO  (usually  shells),  the 
response  is  significantly  different  at  intermediate  to  late  times  compared  to  solids  of  the 
same  materials.  This  code  has  also  shown  some  of  the  inherent  properties  of  conductive 
target  responses  that  may  influence  fundamental  system  design.  For  example  for  a 
conductive,  permeable  steel  sphere,  the  transient  responses  for  the  field,  B,  and  the  time 
derivative  of  the  field,  dB/dt,  are  plotted  in  Figure  A.  1.4.1  below.  Over  six  decades  of  time 
B  falls  by  only  3  orders  of  magnitude,  whereas  dB/dt  falls  by  over  6  orders  of  magnitude. 


Figure  A.l.4.1:  Transient  responses  for  B  and  dB/dt  for  6  cm  conductive,  permeable  steel 
sphere,  at  1  m  depth. 

The  other  important  result  seen  in  Figure  A.  1.4.2,  for  the  same  sphere  but  for 
varying  shell  thickness,  is  that  the  dB/dt  response  becomes  exponential  at  about  10  msec 
for  a  thin  shell,  whereas  it  becomes  exponential  almost  a  decade  earlier  in  B.  From  an 
identification  point  of  view,  the  properties  of  the  shell  are  identifiable  at  earlier  times  and  at 
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higher  signal  levels  in  B  than  in  dB/dt.  Identifying  a  target  through  its  spectral  response 
clearly  must  be  done  with  a  shell.  To  our  knowledge,  no  programs  currently  exist  that  can 
model  spheroidal  or  other  shells. 


Figure  A.l.4.2:  Transient  responses  for  B  and  dB/dt  for  6  cm  conductive  permeable  steel 
sphere,  for  varying  shell  thickness,  at  1  m  depth. 


We  have  built  a  graphical  user  interface  for  this  modeling  code  to  make  it  more 
accessible  to  the  UXO  community.  A  snapshot  of  the  interface  is  shown  in  Figure  A.  1.4.3. 
(Parties  interested  in  using  this  software  need  to  send  a  formal  request  to  LBNL.)  The 
graphical  user  interface  allows  the  user  to  model  a  response  of  either  a  solid  sphere  or  a 
spherical  shell  at  a  specified  depth.  The  source  can  be  a  dipole,  a  plane  wave,  or  a  loop. 
The  receivers  are  vertical  dipoles  that  can  be  placed  in  any  arbitrary  position.  The  response 
can  be  calculated  either  in  frequency  or  time  domain.  For  a  stand-alone  configuration,  the 
program  will  calculate  a  response  over  a  specified  range  of  frequencies  or  times.  In  a 
profile  configuration,  a  response  is  calculated  at  one  particular  frequency  or  time  along  the 
length  of  the  profile.  Time  domain  calculations  allow  for  a  “step  off,”  “step  on,”  or  “square 
wave”  transmitter  current  waveform  signal.  Moreover,  the  code  has  the  option  to  calculate 
a  dB/dt  response  as  seen  through  an  inductive  receiver  of  a  specified  center  frequency  and 
bandwidth. 
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Figure  A.l.4.3:  Graphical  user  interface  for  calculating  responses  of  a  spherical  target 
with  an  arbitrary  transmitter-receiver  configuration. 

A.1.5  3-D  half-space  of  arbitrary  conductivity  and  permeability 

We  have  obtained  the  general  finite  difference  3D  code  developed  by  Greg 
Newman  at  Sandia  National  Labs  (Newman  and  Alumbaugh,  1997).  We  will  use  this  code 
to  model  the  effects  of  inhomogeneous  ground,  particularly  inhomogeneities  in  magnetic 
susceptibility,  on  loop-loop  systems.  This  code  was  recently  modified  by  Newman  and 
promises  to  greatly  accelerate  our  analysis  and  understanding  of  the  role  of  geologic  noise 
in  AEM  systems  for  UXO. 
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